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Brief  Background 


Per  the  description  of  the  Statement  of  Work  of  GRANT  No.  N00014-95-1-G037.  this  two-year 
grant  was  intended  to  fill  a  number  of  current  and  near-future  needs  of  NRL's  cosmic-ray  basic 
research  efforts  as  well  as  related  applied  research  effort?  like  SEE  (space-radiation  environincni 
studies)  that  have  been  developed  or  currently  being  extended  by  the  cosmic  rays  group  at  NR  L. 

The  support  has  been  analytic  and  computational  in  nature  in  the  general  area  of  modeling  of 
cosmic-ray  transport  in  the  heliosphere. 

This  final  report  describes  the  main  support  afforded  by  this  grant  to  NRL’s  cosmic-ray  group 
for  the  two  years  (9/95-9/97)  of  the  grant.  The  main  task,  by  far.  has  been  the  development,  testing, 
and  data  comparison,  of  a  global  time-dependent  and  three-dimensional  heliospheric  transport  code 
of  galactic  cosmic  rays.  While  the  code  is  based  on  current,  standard  and  established  theory  of 
solar  modulation  of  galactic  cosmic  rays,  it  is  far  more  computationally  demanding  thaji  what  a 
typical  application-oriented  .study  may  require  (e.g..  tasks  with  connection  to  SEE  work).  To  such 
ends,  the  purpose  of  developing  such  a  code  wa.s  to  afford  the  group  a.  fully  three-dimensional  solar- 
modulation  model  by  which  a  computationally  efficieni  parametric  set  of  simulated  data  (that  can 
easily  and  efficiently  be  incorporated  in  larger  semi-empirically  based  models  like  CREME96  (Tylka 
et  al.  1997)  for  example)  can  reliably  and  efficiently  be  developed  (e.g.,  .4dams  ic  Lee  1996). 

This  development  should  afford  the  group  lx)th  qualitative  and  quantitative  advantage  in  this 
regard  and  relative  to  modulation  codes  currently  available  to  the  group  which  tend  to  be  rutli- 
mentary  one-dimensional  (so-called  spherically-symmetric)  codes. 

The  first  section  of  this  report  describe, s  in  some  detail  the  basic  physical  model  the  new  three- 
dimensional  transport  code  is  ba.se<l  upon,  'fhe  second  section  highlights  the  numerical  imple¬ 
mentation  and  various  algorithms  of  the  code,  while  t  he  main  routines  are  listed  as  an  appeiidi.'c. 
The  third  section  illustrates  .some  sample  calculations  and  comparison  to  available  data.  The 
fourth  section  suggests  some  directions  and  recommendations  for  future  work  ba.sed  on  this  code 
for  ai>plications-oriented  studies  by  the  group  at  NRL.  as  well  as  for  more  basic-physics  oriented 
improvements.  The  fifth  and  final  section  is  a  list  of  citerl  references. 
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1.  Solar  Modulation  of  Galactic  Cosmic  Rays 

1.1  The  Physical  Model 

The  standard  model  of  long-term  solar  modulation  of  GCR  that  is  iise<l  for  our  purposes  in  iliis 
work  is  the  one  based  on  Jokipii  Ac  Parker  (1970)  transport  e(|uatioii: 

^  =  V  .  (K>VU)  -  V  •  -  (V^) .  rf-  -P  i(V  - A(ar(r)  .  (1) 

where  U  =  U(7\T,t)  is  the  number  density  of  GCR  as  a  function  of  position  r,  kinetic  energj’  T, 
and  time  t.  U  is  related  to  the  (solar-minimum)  obsermble  omni-directional  cosmic-ray  intensity  j 
ns  j  =  rf7/(4r).  where  v  is  the  particle's  speed.  The  streaming  flux  vector  F  is  then  written 

F==-K^Vr  +  {V,^r  +  V,^:-~(aTU)]  .  (2) 

and  the  anisotropy  vector  is  ^  =  F/vU.  Below  we  briefly  describe  the  various  terms  in  E<<.  (I) 
and  the  associated  physical  processes  they  represent,  o  in  Eqs.  (1)  and  (2)  is  tlie  standard  scalar 
function  a  =  (2moC'  -t-  T)/{mo(r  +  T),  with  iiiaC'  being  the  rest-mass  of  the  particle. 

Eq.  (1)  is  a  dynamic-balance  (i.e.,  time-dependent)  statement  for  U  in  three  dimensions  that 
is  subject  to  the  four  fundamental  physical  processes  that  comprise  the  standard  model  of  so¬ 
lar  modulation:  (l)diffiision  (due  to  the  irregular  component  of  the  heliospheric  magnetic  field): 
{2) convection  (due  to  the  outflowing  solar-wind  plasma  carrying  with  it  the  frozen-in  heliospheric 
magnetic  field  lines):  (Z)dnft  (due  to  the  large-scale  curTOture  and  gradient  of  the  regular  compo- 
iieiii  of  the  heliospheric  magnetic  field,  and  {4)ndial/ottc  energy  loss  (due  to  the  diverging  solar-wind 
plasma). 

The  heliospheric  magnetic  field  is  taken  to  be  composed  of  an  irregular  component  superposed 
on  a  regular  one.  i.e.,  B  =  Bo  +  HB.  with  V  •  S  =  0.  We  use  the  standard^  description  for  the 
regular  component  of  the  heliospheric  magnetic  field,  i.e..  a  multi-sector  field  with  Parker’s  spirals 
aiul  with  opposing  polarities  above  and  below  a  wavy  current  sheet  (Kota  Ac  .lokipii  1983:  and 
references  therein): 


Boir.  e.  o)  =  B,Ar.  d)  [J  -  2S{d  -  d')]  .  (3) 

where  the  single-polarity  field  Bo{r.d).  note  its  (viiidepeiidence.  is  written 

Bo(i%  d)  =  .4  [1  »•  -  n,;,  sin  d/\A,,r  .  (4) 

*  .As  ol'  yel-.  our  iiuinerical  code  does  not  take  advantage  of  a  recent  suggestion  (Jokipii  fc  Kota  1989;  Smith  Biel>er 
1991)  ;is  well  its  observation  (Smith  et  cil.  199.5)  to  modify  the  average  field  in  the  polar  regioits.  This  modification, 
however,  can  easily  l>e  made  as  well  as  others  its  outlined  in  ‘j-1.2  of  this  i-eport. 
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(/•,  0, 0)  refers  to  a  heliocentrit  spherical  coordinate  system  that  we  adopt  here.  .4  is  a  constant  that 
carries  both  the  strength  and  polarity  (depending  on  the  solai-  cycle)  of  the  single-polarity  field, 
and  f2<j,  is  the  rotational  speed  of  the  Snn  around  its  axis.  V',„.  is  the  solar-wind  speed  with  a  radial 
but  ^-dependent  velocity  profile*  ; 

=  v;(l  -I-  /^sin'  X)r  ,  (5) 

where  A  is  the  solar  geomagnetic  latitude.  is  the  solar-wind  speed  at  zero  latitude,  and  //  is 
;i  fitting  parameter  correlated  with  the  solar  cycle  (.Ananthakrislman  et  al.  199-5).  .Note  that 
V  •  I  7^:  0.  S(0  —  0')  is  a  step  function  with  0'  given  by 

—  -I-  .siii“‘  [sin  n  sin(c>  —  o,-,  -f  rn.j,/V’su,)]  .  (fi) 

and  a  is  another  fitting  parameter  related  to  the  tilt  angle  of  the  current  sheet  at  the  Sun  (also 
correlated  with  the  solar  cycle),  and  6,-,  is  an  arbitrary  constant.  The  current  sheet  is  essentially 

defined  by  S(^-«')-^ 

The  drift  velocity  vector  averaged  over  near-isotropic  pitch-angle  scattering  is  written  (Jokipii 
ot  al.  1976): 


=  (7) 

where  P  is  the  magnitude  of  the  particles  momentum  with  charge  q  and  i3  =  t'/c,  where  c  is  the 
speed  of  light.  Note  that  V  •  (Vj)  =  0,  and  {Vj)  is  both  polarity  and  charge  sensitive  via  qA  being 
>  0  or  <  0.  Angle  brackets  in  (i  rf)  denote  averaging  of  pitch-angle  resonant  scattering  along  the 
field  lines  due  to  the  irregular  component  of  the  field. 

'Phe  irregular  component  of  the  magnetic  field  is  taken  to  be  a  random  stationary  field  with  zero 
mean.  {SB)  =  0.  and  characterized  by  a  Kolmogorov-like  power  spectrum  (i.e..  oc  !F{{8B  •  oc 
/.•“*•)  (.lokipii  1971).  with  k  being  the  wavenumber.  C,  the  Kolmogorov  spectral  index,  and  T{-  ■  •} 
denotes  Fourier  transform.  The  symmetric  diffusion  teu.sor  in  the  local  solar-wind  frame  is  written: 

/ 0  0  \ 

A’*  =  0  Kx  0  .  (8) 

Vo  0  Kii/ 

■  Unlike  tlie  discussion  in  Note  1.  and  althou)i;h  this  dependence  was  not  typically  part  of  the  standard  modulation 
model,  ii  was  included  in  our  numerical  code  eaily  on  due  to  the  full  time-dependent  three-iliniensional  physical 
pin  lire  it  was  built  to  model.  A  second  reason  for  its  eiU'ly  iiictu.sion  had  to  do  with  our  paying  special  attention 
to  the  question  of  inner  boundaiv  contlilions  and  the  so-calleil  Peclet-sjieclruin  analysis  that  we  perform  in  oiir 
numerical  solution  of  the  convective-diffusive  Eq.  (1).  *>2.2 

■^See  accomptinying  figure  for  a  perspective  of  the  magnetic  cun-ent  slieet  in  cartesian  coordinates. 
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where,  from  quasi-linear  theo^\^ 


K||  =  Ko  R‘  ^  3 


(9) 


with  E  being  the  particle's  rigidity.  is  the  strength  of  the  field  at  Earth,  and  Ho  is  a.  constant.  The 

♦ 

la.st  term  in  Eq.  (9)  arises  from  recjuiring  that  the  particle’s  gyroradius  remains  .smaller  than  t  lie 
•scatlering  mean-free- path  (Morfill  tV.  Volk  1979;  .lokipii  ic  Davila  1981)  throughout  the  modulation 
region.  We  will  have  more  to  say  about  kx  in  ‘jl.'i.  In  the  heliocentric  si)herical  coordinate  sysietn 
A'*  becomes 


where,  in  terms  of  K||  and  kx. 


Hyy  0  \ 

0  0  I 

^0r  .0  -  j 


K,.,  =  /v||  COS*  -i-  Kx  sin*  . 
K,,,.,  =  (kx  -  K||)cos4'sin  . 
'>■■«(>  =  • 

K/fty  —  ^-^0  * 

««>e  =  «x  cos*  ^  -b  K.||  sin*  ^  . 


(10) 


(11) 


a  nd  ']'  =  tan  ’  (rfi,t./V’,u.). 

I'o  complete  the  physical  picture  of  the  standard  modulation  model  one  needs  to  specify  phys¬ 
ically  meaningful  initial  and  boundary  conditions.  We  will,  however,  defer  this  discussion  ti>  h'..’.-- 
where,  due  to  the  full  time-dependent  and  three-dimensional  nature  of  the  model  and  its  numericai 
implementation,  coupled  with  a  rather  stiff  and  complex  transport  PDE.  we  pay  special  attention 
to  such  conditions  in  our  Peclet  analysis  section  of  the  transport  process. 


1.2  Standard  Queisi-Linear  Cross-Field  Diffusion 

The  theory  of  cros.s-field  diffusion  and  its  applications,  e.g..  to  cosmic  rays  heliospheric  transport, 
owes  its  inception  to  the  works  of  .lokipii  and  Parker  (.lokipii  1966:  .lokipii  and  Parker  1969a;  .lokipii 
197.3)  where  the  mechanism  was  understood  as  mainly  a  non-resonant  one  due  to  the  random-walk 
of  tile  magnetic  field  lines  themselves,  so  that  t  est  particles  propagat  ing  along  the  lic'ld  lines  will  also 
diffuse  normal  to  the  lines.  Cross-field  diffusion  due  to  resonant  .scattering  was  shown  to  contribute 
little  to  the  mechanism  (Jokipii  k  Parker  1969a.).  [.Note  that  the  same  mechanism  was  also  being 
studied  in  the  context  of  plasma  physics  and  applications  to  Tokamak-field  turbulence  (Ro.senbluth 
et  al.  1966:  for  a  recent  review  see.  e.g..  Isichensko  1992).] 
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In  addition  to  application  to  cosmic  rays  heliospheric  transport,  which  is  the  focus  of  this  work, 
cross-field  diflfusion  also  plays  significant  roles  in  other  astrophysical  applications  like  diffusive  shock- 
acceleration  theories  (Drury  1983;  Lagage  &  Cesarsky  1983;  Blandford  &  Eichler  1987:  .lokipii  1987: 
.lokipii  fc  Morfill  1987:  Achterberg  ic  Ball  1994;  Duffy  et  al.  1995:  Giacalone  k  .Jokipii  1996)  and 
cosmic  rays  galactic  transport  (Jokipii  k  Parker  196§b;  Skilling  et  al.  1974:  Ptuskin  1979:  Barge 
et  al.  1984:  Corbelli  &  Veltri  1989:  Chnvilgin  Ptuskin  1993:  .Achterberg  A:  Ball  1994:  Cliacalone 
Jokipii  1994:  Dendy  et  al.  1995;  Duffy  et  al.  1995;  Klepach  1995;  Ptuskin  1995). 

The  two  essential  ingredients  in  the  theory  of  cross-field  diffusion  are  the  separation  rate  of  the 
two  (initially  close)  wandering  field  lines  and  the  scaling  of  the  perpendicular  diffusion  coefficient 
witli  the  strength  of  the  magnetic  turbulence.  The  separation  rate  is  understood  to  be  a  critical 
length  scale  in  the  problem  as  the  two  lines  (and  test  particles  that  are  tied  to  them)  tend  to  become 
independent  beyond  this  length  scale.  Early  estimates  of  the  separation  rate  (.lokipii  1973:  Ptuskin 
1979:  Barge  et  al.  1984)  give  a  rate  on  the  order  of  the  inverse  of  the  turbulence  <'orrelalion  length. 
[More  recent  numerical  (Zimbardoet  al.  1995:  Zimbardo  fc  Veltri  1995)  and  analytic  (Barghouty  A 
Jokipii  1996)  estimates  suggest  a  more  complex  separation  rate  that,  in  addition  to  the  correlation 
length,  is  coupled  to  the  strength  of  the  turbulence  in  a  non-trivial.  multi-region  fashion.  See  §4.2 
for  more  on  the  role  of  nonlinear  cross-field  diffusion,  as  opposed  to  quasi-linear  description,  in 
GC'R  modulation  modeling.] 

The  standard  quasi-linear  description  (Jokipii  1971)  for  the  random  walk  of  the  field  lines,  in 
the  limit  when  the  scattering  mean  free  path  is  ^  the  turbulence  correlation  length  (..  is  expre.ssed 
as 


U-r-) 

A- 


tw(0,0.r') 


VrAI'-  =  0) 


(12) 


with  Crr  being  the  correlation  function  of  SBj.  taken  in  the  limit  Ar  »  L  and  T’xx(A:  =  0)  is  the 
corresponding  power  spectrum  at  zero  wavenumber,  kx  Is  then  given  by 


S'x  ~  ■j^'PrAI'-  =  0)  .  (13) 

[The  (|uasi-linea.r  scaling  is  evident  here  in  that  the  diffusion  coefficient  scales  with  Sir.]  Comparing 
t  he  above  expression  to  a  similarly  derived  expres.sion  for  K||  al  high  rigidities 


K.|| 


rs./ 


VxA^'  —  0) 


(14) 


and  using  an  extrapolated  estimate  for  VxA^  =  0)-  is  led  to  qualitatively  conclude  that 


Kx  ^'ii  • 


(lo) 


Final  Rcpori /NOOOl  1-!).V l-G();ir/r> 


\;  ept.  perhaps,  for  very  small  rigidities  where  the  al>ove  expressions  do  not  apply. 

Thus  in  standard  modulation  models  of  GCR  it  has  become  a  standard  practice  to  assume  that 
the  magnitude  of  kx  is  simply  only  a  small  fraction  of  K||.  typically  ^  10%,  i.e.,  and  apart  from  the 
actual  numerical  ratio  one  would  use  for  /vx/''(|-  '>x  's  essentially  couplcKl  to  K||.  Therefore,  and 
while  we  do  not  address  more  recent  suggestions  for  variants  of  this  assumption  (e.g..  anisotropy  in 
Kx,  modifications  due  to  larger  fields  at  high  latitudes,  etc.  (see  discussions  in.  e.g..  Kota  &  Jokipii 
1995.  Bieber  et  al.  1995),  for  our  purposes  here,  it  must  be  kept  in  mind  that  the  quasi-linear  scaling 
is  wliat  is  typically  assumed  in  standard  modulation  models,  as  in  the  above  standard  quasi-linear 
expressions,  as  well  as  in  variations  thereof,  independent  of  the  assumed  numerical  ratio  sx/''i|- 
In  arldition  to  the  standard  quasi-linear  scaling  alluded  to  above,  the  staiulard  .separation  rate 
as  deduced  from  quasi-linear  theory  (.lokipii  1973:  Barghoutv  tk:  .Jokipii  1996)  is  dependent  l^otli 
on  the  so-called  large-step  diffusion  coefficient.  =  (A.r-)/Az  as  given  in  Eq.  (12)  above,  with 
its  quasi-linear  scaling  embedded  in.  as  well  as  on  the  parallel  length  scale  in  t  he  turbulence  T  To 
illustrate  this  for  a  turbulent  field  with  truncated  Kolmogorov-like  power  spectrum,  i.e., 


P i;Tik) 


t)=  (l^vex 


exp{ik  •  c^^.(f)  a 


(kl  +  kl) 

(A-2  -I-  f-2)C/2+-> 


exp(-fr'ff)  . 


(16) 


with  /.•  being  the  wavevector  and  C;  (•<  f)  being  the  turbulence  inner  scale,  the  Fokker-Planck 
diffusion  coefficient  Drp(p)  =  "here  p  =  |^rx|  C  (  and  (Jrx  =  (Sx.Sy)  is  the  2-d 

separation  vector  between  the  two  lines  in  the  r-y  plane,  can  —to  first-order  in  =  {(■,/()’—  b<' 
written  as  (Barghouty  Jokipii  1996) 


Dfp(p)  «  2£>oc. 


1  - 


2 


r(a2)  in) 


c/2 


(17) 


where  1\„{Z)  is  the  Bessel  function  of  an  imaginary  argument  and  F  is  Euler’s  gamma  function,  t] 
and  t.>  are  expansion  coefficients  in  d('  given  l>y 


=  1  -h 


<U- 

(2-(  +  -isr^}  ' 


<>f-’(3  +  C) 

(2-(  +  -W'f-)  ■ 


In  short,  by  standard  parameterization  of  the  diffusion  coefficient  and  the  separation  rate  in  the 
qua.si-linear  description  of  the  random  walk  of  the  field  lines  we  are  referring  to  a  scaling  oc  Sb- 
and  a  separation  rate  that  depends  on  a  parallel  correlation  length  of  SB.^  Further,  the  statistics 
is  inherently  Gaussian  and  the  geometry  is  slab-like.  Taken  together,  the  quasi-linear  description 


'  Fo]-  a  magnetic  field  in  X.y.  Z  cai'tesian  coordinates  with  a  imiforni  part  B,-,  along  Z  and  a  normal  fluctuating  one 
SB{x.  y.  z).  B  =  BoZ  -h  SB.  with  (SB)  =  0  aii<l  suhjeci  to  V  •  B  =  0.  the  relative  turbulence  strengtii  is 
tlefined  as  Sb  —  {{SB  •  SB})^^~ /Bo-  where  angle  brackets  denote  ensemble  averaging. 
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of  cross-field  diffusion  is,  therefore,  inherently  a  classical  brownian-motion  description  of  the  field 
lines. 

2.  Numerical  Implementation 

2.1  Numerical  Solution  of  the  Transport  El^uation  in  l-f4  Variables 

We  have  develoi>ed  our  own  numerical  solution  (Barghouty  1997)  that  integrates  over  t  Eq.  (1) 
in  three  spatial  dimensions  (r,9.d>)  and  the  kinetic-energy  variable  T.  i.e..  a  full  numerical  solut  ion 
in  1-f-l  variables.  In  this  section  we  briefly  describe  the  method  by  which  our  numerical  solutio:i  i- 
arrived  at  since  it  deviates  from  some  of  the  earlier  developerl  ones  for  this  problem  (we  are  referring 
here  to  solutions  based,  e.g..  on  Crank- .N'ickolson  techni(|ues  and  alternating  direction  implicit  (or 
.“VDl)  schemes  with  either  momentum  or  energy  as  pseudo-time  (e.g..  Fisk  1971.  1976:  Perko 
Fisk  1983).) 

Our  numerical  code  integrates  Exp  (1)  using  the  so-called  numerical  method  of  lines  (Scheisser 
1991)  in  which  the  time  variable  is  kept  "continuous"  and  the  finite-differencing  scheme  is  applied 
to  the  other  four  variables  at  different  points  along  a  “time-line’*  which  is  then  integrated.  The 
advantage  from  separating  the  time  variable  from  the  rest  is  that  the  stability  issue  of  the  solution 
is.  from  the  outset,  separated  from  the  accuracy  issue.  For  a  PDE  like  Eq.  (1).  simple  eigenvalue 
tests  can  reveal  the  stiffness  of  the  eqtiation  (this  stifl'ness  can  also  be  readily  appreciated  by  t  h<' 
tiiany  orders-of-magnitude  separation  in  the  strengths  of  the  various  diff  usion-tensor  terms  at  widely 
separated  radial  points),  so  that  the  stability  issue  is  of  concern. 

.411  first-order  and  second-order  diffusion-related  derivatives  of  the  dependent  variable  U  on  the 
UH.S  of  Eq.  (1)  are  evaluated  at  a  certain  t  using  a  five-point  centered  differencing  scheme  (i.e..  ac¬ 
curate  to  fourth-order  in  the  .step  size  of  the  respective  variable).  .411  first-order  convection-related 
derivatives  are  evaluated  using  a  five-point  either  upwind  or  downwind  biased  (depending  on  the 
variable)  differencing  scheme.  .411  field,  diffusion  tensor,  and  drift  vector  terms  along  with  their 
first  and  second  order  derivatives  arc'  evaluated  analytically.  Once  all  derivat  ives  of  /'  all  collected 
as  the  RHS  of  E<|.  (1)  at  all  t  points  along  the  time-line,  they  are  integrated  using  the  (truly  un¬ 
conditionally  stable  and  extremely  efficient)  sparse  .Jacobian-matrix  technique  of  Hindmarsh  (1982) 
(see  also  Byrne  Hindmarsh  1987)  appropriate  for  stiff  PDEs.  The  accuracy  in  this  integration 
scheme  is  completely  under  the  control  of  the  user. 

2.2  Boundary  Conditions  and  Peclet  Analysis 

"initial  condition"  for  the  purpose  of  reaching  a  stationary  .solution  means  at  /  =  0  we  assume' 
that  ['  takes  on  the  unmodulated  LISM  isotropic  distribution  at  the  a.ssum<'d  outer  heliospheric 
boundary  of  =  100  Alb  of  the  the  form  ,  0.  o.  T)  =  =  const  -  +T]~'''  where  ;/  is 

the  typical  LISM-spectrum  exponent  2.65-2.75.  This  "initial  condition"  a.ssumption  also  serves  as 
the  boundary  condition  for  U  at  r^.  =  100  .4U.  The  inner  heliospheric  boundary  condition  at  i\,  = 
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.01  AU  is  either  taken  to  be  U{ro,  9. T)  =  0,  or  the  number-density  conserving  condition  according 
to  Liouville’s  theorem,  i.e.,  F(ro,0,0.T)  =  0  (which  was  assumed  in  the  sample  calculations  of  !{.3) 
depending  on  the  level  of  analysis.  For  Eq.  (1)  this  condition  is  expressed  as 


du 

dt 


v; 


(19) 


Both  conditions  seem  to  give  very  similar  results  at  1  AU,  but  this  is  not  a,ssure<l  to  be  the 
case  throughout  the  modulation  region  nor  for  all  energies.  For  boundary  conditions  in  the  angle 
variables,  dUfdO  =  0  =  dUfdo  are  a.ssumed  periodic  over  the  range  [0,  r]  for  9  and  [0,27r]  for  o. 
Because  of  the  complete  symmetry  of  the  .solution  aixtut  the  current  sheet,  the  range  for  9  is  taken 
to  be  [O.r/2]  thereby  treating  the  inner  singularity  in  the  solution  (when  B  -»■  0  at  all  points  on 
the  sheet)  as  a  boundary  condition  for  9.  The  singularity  at  ^  =  0“^  is  treated  using  THospital  rule. 
For  the  “boundary"  condition  in  T  we  simply  assume  that  U(r.9,<i).Ty^.)  =  U(){r,^,.9.o.  T^,).  with 
7'x.  =  100  GeV'/Nucleon. 

We  perform  a  Peclet  analysis  of  Eq.  (1),  and  as  we  briefly  describe  below,  t  he  analysis  touches 
upon  both  the  viability  of  the  numerical  solution,  on  the  one  hand,  as  well  as  the  assumed  boundary- 
conditions.  on  the  other.  Peclet  number  (or  spectrum  for  energy-dependent  transport  coefficients 
as  in  Ek].  (1))  is  defined  as  the  ratio  of  the  convection  group  to  the  diffusion  group  in  the  PDE 
times  the  characteristic  length  scale  in  the  problem.  For  Eq.  (1),  for  example,  w-e  can  define  the 
Pwlet  spectrum  for  the  9  and  r  convection  and  diffusion  groups  as 


V,{9,T)  =  |(i‘^)«(r”.#,d-.r)|-/WK,«(/--.6>.o-.r)  .  (20) 

VA9.  T]  =  |(K/),.(r“.  o'.  T)  -F  V;„.(6>)|r,^,/K,. 9.  o'.  T)  .  (21) 


for  some  /■'  and  o'-  The  significance  of  this  number  is  that  when  it  is  too  high  1),  indicativ(' 
of  a  strongly  convective  PDE,  one  needs  to  pay  special  attention  to  the  behavior  of  the  numerical 
solution  at  both  exterior  and  interior  boundaries  as  it  can  admit  spurious  discontinuities  which 
can  propagate  throughout  the  characteristic  length.  In  other  words,  the  PDE  can  (numerically)  be 
iiuide  to  resemble  more  an  advection  ecpiation  with  its  usual  reflective  properties,  but  with  spurious 
implications  for  its  diffusive  properties. 

In  numerical  studies  of  Eq.  (1).  and  especially  with  issues  related  to  diffusion  in  three  dimen¬ 
sions.  this,  in  turn,  may  lead  to  erroneous  and  numerically-induced  spurious  “pl'.vsicf^-l  features  in 
the  transport  (e.g.,  Schiesser  1996).  To  illustrate,  in  Fig.  I  we  plot  Vs(9.T)  for  H  calculated  at 
1  AU  but  averaged  over  p.  What  Fig.  1  seems  to  indicate  is  that  except  for  high-energy  protons 
and  close  to  the  sheet  in  the  ecliptic,  the  ^-transport  is  weakly  convective  (or  mostly  diffusive). 
C'ontrast  this  with  the  'Pr{9,T)  spectrum  in  Fig.  2  for  the  radial  groups  where  for  low-energy 
protons  and  essentially  for  all  9  the  7-transport  is  primarily  convective. 
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3.  Sample  Calculations  and  Comparison  to  Data 

To  give  a  feel  for  the  high  efficiency  of  the  integration  algorithm  we  use.  for  the  sample  calcu¬ 
lations  we  show  below  and  using  a  grid  size  of  lo  x  15  x  15  x  15  (radial  and  energy  grid  points  are 
logarithmically  spaced),  i.e.,  50.625  ODEs  to  be  solved  in  general,  and  with  a.  prescribed  relative 
accuracy  of  1%  (for  time-integration),  the  method  required  the  actual  evaluation  of  only  11.99S 
ODEs  (or  grid-point  visits)  due  to  the  sparsity  of  the  .lacobian  matrix  of  the  ODE  system  (about 
24%  of  total  number  of  ODEs  that  other  methods,  generally  speaking,  may  liave  to  solve).  With 
increasing  grid  size,  the  efficiency  gets  even  better  as  the  sparsity  increases  and  an  order  of  mag¬ 
nitude  increase  in  efficiency  is  not  untypical.  Reaching  a  stationary  solution  required  a  mere  67 
points  along  the  time-line  for  a  total  CPU-elai)sed  time  of  3.1  hrs  (or  about  23  ODEs  per  sec)  on 
a  DEC-ALPHA  250  machine. 

For  illustration  purposes  (our  focus  is  on  .solar-minimum  .steady-state  solniions)  Ecp  (1)  is 
integrated  until  a  stationary  solution  is  reached  (/  ^  0.15  yr)  with  the  prescribed  a.c<'nra.ry.  Tinic- 
sensitive  solutions  are,  naturally,  also  available.  On  a  relatively  coarse  grid.  Fig.  3  shows  a  satnple 
calculated  (i-averaged  H  flux  at  1  AU  and  B  =  90^  (in  units  of  particles/m --sr-s-MeV/Nucleon) 
at  solar  minimum.  Data  from  the  1977  and  1987  solar  minima  are  also  shown  (even  though  the 
calculations  pertain  to  >  0  ,i.e..  1977  solar-minimum).  Table  1  lists  the  salient  assumed  physical 
parametei-s  in  the  illustrations.  Fig.  4  shows  the  calculated  solar-minimum  H  flux  as  a  function  of 
the  polar  angle  B  at  1  AU  for  three  different  energies  where  relaxation  to  isotropy  with  increa.sing 
energy  is  evident.  In  Figs.  5  and  6  the  local  B  and  r  gradients  are  plotted  for  H  al  solar  minimum. 

4.  Recommendations  for  Future  Work 

4.1  Semi-Empirical  Modeling 

Semi-empirical  expressions  for  an  effective  long-term  modulation  of  galactic  cosmic  rays  can  be 
put  forward  under  the  assumptions  of  steady-state  and  spherical  symmetry  conditions  as  has  been 
done  by.  e.g..  Evenson  et  al.  (1983).  Carcia-.Munoz  et  al,  (1985),  and  more  recently  by  Adams  A 
Lee  (1996)."  In  these  semi-ein])irical  expressions  convection,  diffusion  and  adiabaiic  energy  losses 
[but  no  drift  effects]®  are  all  lumped  into  an  effective  diffusion  coefficient  /!(/•.  T).  .so  that  oim'  •  ;iii 

■’One  can  even  use  an  eflFective  modulating  potential  (ttleeson  .y  .Axford  1968)  in  a  numerical  siniulation  as  has  been 
clone  by.  e.g..  Letaw  et.  al.  (1984). 

"It  is  possible  to  heuristiccilly  include  sucli  effects  in  this  seiui-einpirical  context  by  lumping  the  radial  component 
of  the  drift  velocity  vector  to  the  radial  solar-wind  velocity  vector.  However,  this  drift  addition  is  due  to  all  spatitd 
inhoinogeneities  in  the  regular  heliospheric  magnetic  field  and  not  just  along  the  radial  direction.  As  such,  and  in  a 
spherically  symmetric  picture  of  the  modulaiion  region  one  is  not  justified  in  . including  only  the  i-adiol  componeni 
of  tlio  cli'ifi  while  ignoring  the  others! 
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arrive  by  direct  integration  at  a  simple  expression  of  the  form 

U{r.T)  =  Uo{T)exp 

where  the  spatial  dependence  of  D{r,T)  is  separated  from  its  T  (or  rigidity)  dependence  and  is 
made  to  account  for  the  adiabatic  energy  losses. 

Semi-empirical  expressions  like  the  one  above  appear  to  be  insensitive  to  the  assumed  outer 
l>oundary  of  the  modulation  region  as  well  as  to  the  exact  form  by  which  D{r.  T)  depends  on  both 
;■  and  T.  .41so.  such  expressions  appear  to  give  reasonable  fit  to  near  Earth  data  over  the  1  l-\  i 
solar-cycle,  i.e..  not  just  for  solar-minimum  conditions,  when  D  is  allowed  to  take  on  different  values 
over  t  he  cycle,  i.e..  taking  D  as  function  of  time  in  addition  to  r  and  T.  .Modulation  studies  using 
such  semi-empirical  expressions  tend  to  converge  to  the  main  conclusion  that  ,  in  solar  a.s  jiear  Earth 
data  are  concerned,  the  steady-state,  spherically-symmetric  modulation  model  of  galactic  cosmic 
rays  is  not  an  unreasonable  one. 

[For  near  Earth  data  one  needs  to  keep  in  mind  that  this  regime  of  the  overall  modulation  region 
represents  but  a  deep  and  hence  asymptotic  (in  the  mathematical  sense  of  the  transport  equation 
and  in  the  sense  of  the  distance  that  the  galactic  c(»mic-ray  particle  has  to  travel  from  the  outer 
boundary  of  the  heliosphere  to  near  Earth)  regime  of  the  solution  both  in  space  and  time.  Thus, 
it  may  not  be  all  that  surprizing  that  the  steady-state  assumption  appears  to  be  a  rea.sonable  one. 
riie  si)herical-symmetry  a.ssumption.  on  the  other  hand,  is  less  obvious.  But  one  can  still  argue 
l  hat  at  lower  energies  (^  GeV/.Nucleon).  where  the  modulation  is  much  more  pronounced  and 
significant,  adiabatic  energy  losses  dominate  over  the  diffusive  (and  other)  aspects  of  the  transport 
|)rocess  and  as  a  result  the  modulation  is  less  sensitive  to  variations  in  the  diffusion  terms  (as  well 
a.s  assumed  geometry),  especiallv  so  for  near  Earth  oijservations.] 

One  possible  semi-empirical  direction  one  can  take  is  to  use  the  simulated  results  of  the  :JL) 
code  to  build  a  data  base  of  the  simulated  modulating  scalar  factor  U(t~.  T)/Uo[T)  throughout  the 
heliosphere  and  as  a  function  of  rigidity.  This  way  any  transported  particle  at  any  given  point  cis 
a.ssigned  a  factor  depending  on  its  rigidity.  Computationally  this  should  lx*  (piite  elticient  .  il  not 
all  that  robust.  Semi-empirically.  one  can  think  of  this  modulating  scalar  funct  ion  as  similar  to  t  he 
one  appearing  in  Eq.  (22)  in  the  lowest  order,  i.e..  in  the  sense  of  expansion  in 

4.2  Improvements  to  the  Physical  Model 

4.2.1  Nonlinear  Cross-Field  Diffusion 

In  recent  studies  of  cross-field  diffusion  (Zimbardo  et  al.  1995:  Zimbardo  .C  Veltri  1995)  and 
analytic  (Barghouty  &  -lokipii  1996)  the  nonlinearity  appears  to  be  an  overriding  characteristic  that 
points  to  the  non-Gaussian  statistics  in  the  problem,  i.e..  the  random-walk  of  t  lu-  held  hues  can 
no  longer  be  viewed  as  a  brownian-like  motion  with  Gaussian-like  statistics  but  rat  her  (dei)endiug 


/ 


drv;„,D-'(r,r) 


(22) 
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on  the  strength  of  the  turbulence)  as  Levy  flights  with  long-range  correlations  according  to  noii- 
Gaussian  (Levy)  distributions  (e.g..  Shlesinger  et  al.  1995).  .4nother  interesting  finding  in  the 
same  studies  is  that  the  statistics  appear  to  be  consistent  with  the  standard  quasi-linear  (whicli 
is  consistent  with  Gaussian)  statistics  when  the  relative  turbulence  strength  is  large  enough.  The 
studies  have  shown  that  for  Sb  ^  .2  quasi-linear  statistics  is  a  valid  description  of  the  random- 
walk  of  the  field  lines,  whereas  for  Sh  ^  .2  the  lines  tend  to  exhibit  non-Ganssian  diffusion  and 
a  c|iiasi-linear  description  appears  inadequate-  Formal  analysis  of  this  behavior  attributes  tlii.s 
iioii-Gau.ssian  diffusion  to  the  effects  of  higher-order  correlations  and  magnetic  percolation. 

Turning  to  the  scaling  of  the  perpend icular  diffusion  coefficieni  with  the  relative  turbulence 
strength,  another  important  facet  in  tlie  description  of  cros.s-field  diffusion,  the  recent  numerical 
calculation  of  Gray  et  al.  (1996)  has  also  shown  that  the  quasi-linear  scaling,  i.e..  D±  oc.  Slr 
appears  adequate  for  large  Sb.  consistent  with  the  non-linear  calculations  alluded  to  above.  In  this 
study  the  fluctuating  field  is  taken  to  be  composed  of  a  slab-like  and  a  2-dimensional  components. 
i.('..  SD{x.i).z)  =  SB-2D(^'-y)  +  <5i5,,(r).  What  was  found  (for  models  appropriate  for  .solar  wind 
1  urijulence)  is  that  the  scaling  D±  oc  Sb  is  the  correct  one  for  an  >tO/ir  2-dimensionaJ  +  20%  slab 
turbulence.  While  the  .study  a,ssuined  classical  diffusion  (i.e..  Gau.ssian  statistics  was  used  in  the' 
numerical  realization  of  the  fluctuating  field),  it  nonetheless  points  to  another  interesting  feat  ure 
of  the  turbulence  in  that  a,s  Sb  —¥  0  Dj_  — >  D>o-  whore  D^o  is  2D-fl net  nations  diffusion 

coefficient  as  described  in  the  non-perturbative  theory  of  field  lines  random-walk  of  .Matthncin^ 
et  al.  (1995).  D^d  —  l/Bo-  where  fx  is  the  correlation  length  of  the  2D  turbulence,  and 

Dx  =  1/2  [Dg  -f  -y/Dj  -f  AD^q).  This  result  is  to  be  contrasted  with  the  slal>-like  coefficient. 

=  SB'l( /'2B'l.  where  C  is  the  correlation  length  along  the  mean  field.  In  short,  it  appears  from 
t  his  study  that  in  the  limit  of  small  turbulence  level,  the  perpendicular  diffusion  coefficient  in  2-d 
t  iirbulonce  geometries  depends  on  a  lengtii  scale  other  t  han  a  parallel  one  (a.s  t  he  C|nasi-iinea.r  theory 
would  suggest),  and  appears  not  to  follow  the  <niasi-linear  .scaling. 

For  the  purpose  of  this  work,  it  appears  from  these  recent  studies  of  the  random  walk  of 
the  held  lines  that:  (/■)  On  the  one  hand.  c|ua.si- linear  description  may  still  be  appropriate  for 
large  relative  turbulence  strength  irrespective  of  the  assumed  turbulence  geometry  (2-d  vs.  slab) 
and  statistics  (non-Gaussian  vs.  Gaussian).  This  finding  is  rather  unexpected  since  in  the  very 
fU'velopment  of  cpiasi-linear  theory  assumptions  about  small  turbulence  levels  were  made  for  the 
theory  to  be  applicable  and  it  was  further  a.ssumed  that  for  large  turbulence  levels  c|uasi-linear 
theory  may  not  be  applicable  due  to  the  method  by  which  the  diffusion  coefficient  was  deriv('d. 
i.e..  in  the  limit  of  Shi-  1  (.lokiini  1966.  1971:  .lokipii  .4.  Parker  i969a).  (/'/)  On  the  other 
hand,  quasi-linear  description  of  field  lines  random  walk  has  been  shown  by  these  studio's  to  Im' 
inadequate  for  small  turbulence  levels  under  the  assumption  of  either  Gaussiiui  or  non-Gaus.sia.ii 
statistics  and  especially  for  2-d  turbulence  geometries.  This  is  a  far  reaching  conclusion,  and  its 
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implications  for  the  heliospheric  transport  of  cosmic  rays  deserves  concerted  analytic  as  well  as 
numerical  efforts.  One  direction  is  to  couple  the  perpendicular  diffusion  coefficient  to  both  the 
relative  turbulence  strength  and  an  appropriate  normal  length  scale,  thus  coupling  geometry  with 
statistics  in  a  description  of  the  random  walk  of  field  lines. 

Tlie  parameterization  of  the  so-called  (effective)  anoqjalous  diffusion  developed  for  magnetic 
t  urbulence  with  strong  anisotropy  (Rosenbluth  et  al.  1966;  Krommes  1978:  Rechester  &  Rosenbluth 
1978:  Kadomstev  k.  Poguste  1979:  Isichenko  1991.  1992)  can  be  used  to  accomplish  such  a  coupling. 
In  this  nonlinear  description  of  the  random  walk  of  the  field  lines,  the  anisotropy  (due  to  K||  »  kx  tts 
in  solar-wind  turbulence  geometries)  has  been  shown  to  significantly  enhance  the  transport  across 
the  mean  field,  i.e.,  giving  rise  to  an  effectively  large  k][. 

For  a  general  2-d  magnetic  turbulence  characterized  by'  correlation  lengths  along  and  across  the 
the  mean  field,  fn  and  fx,  a  Kolmogorov  length  tj^  for  the  exponential  divergence  of  the  field  lines 
in  the  plane  normal  to  the  mean  field,  i.e..  (/>(.:))  =  /),-,e.\p(r/fA-).  where  p(r)  is  as  describerl  in 
Sect.  1.3  and  po  =  p{z  =  0).  and  field-lines  diffusion  coefficient  D,„  (e.g..  D±  in  a  2d-turbuleJic.(' 
geometry  as  alluded  to  in  Sect.  1.1).  the  effective  cross-field  diffusion  is  written 


While  the  above  formulation  was  developed  and  applied  to  magnetic  turbulence  characterized  by 
large  Kubo  numbers,  Tv  =  <56  {(\\/(±)  ^  1.  i-c..  the  so-called  magnetic  percolation  regime,  its 
applicability  does  not  appear  to  be  restricted  to  this  regime  (Isichenko  1991).  We  will  be  applying 
it  in  the  regime  that  is  characteristic  of  space  and  a.strophysical  plasmas,  i.e..  Tv  ~  0{Sb)  ^  1 
(Ziiiibardo  et  al.  1995:  Zimbardo  .V  Veltri  1995:  Barghouty  iV:  .Jokipii  1996). 

For  («;||,  (||)  and  (kx.  (±)  the  paraiiielerization  and  values  of  the  (unenhanced)  standard  descri|> 
tion  are  retained,  while  noting  that  Ky  >  and  (x  —  Tvt||.  For  Z?„,  we  will  take  advantage  of  the 
parameterization  of  Matthaeus  et  al.  (1995).  and  for  C^-  w'e  will  use  a  fit  to  the  entropy  calculation 
of  Barghouty  &  .Jokipii  (1996)  wherein  ( a  —  /2. 

.Note  that  for  small  <56  the  above  parameterization  of  nonlinear  cross-field  diffusion  takes  into 
account  both  the  effects  of  the  2d-t.urbnlence  geometry  on  thediffusivity  of  the  field  lines  (.Ma  tthneii:- 
et  al.  1995:  Gray  et  al.  1996)  via  the  scaling  character  of  £>,„  as  well  as  the  non-Ganssian  .statistics 
of  t  he  lines  random-tvalk  (Zimbardo  et  al.  1995:  Zimbardo  &  Veltri  1995:  Barglionty  <V.  .Jokipii 
1996)  via  l/c-  Moreover,  when  Sb  is  large  enough  the  parameterization  reduces  <|uite  smoothly  to 
the  qua.si-linear  description. 

One  can  first  explore  the  degree  to  which  this  nonlinear  parameterization  deviates  from  the 
standard  one  insofar  as  the  calculated  long-term  cosmic  rays  fluxes  are  concerned  throughout  the 
modulation  region,  with  particular  emphasis  on  latitudinal  transport  and  the  observed  high  degrei' 
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of  isotropy  therein,  e.g.,  Ulysses’  recent  observations  (Mckibben  et  al.  1995a).  Also,  by  calculating 
time-sensitive  cosmic  rays  fluxes,  one  is  able  to  focus  on  short- terms  effects  of  cross-field  diffusion, 
e.g..  interactions  with  CIRs  and  the  current  sheet  (.lokipii  fc  Kota  1995:  Mckibben  et  al.  1995b). 
Finally,  and  from  fits  to  various  observed  radial  and  latitudinal  cosmic  rays  intensity  gradients  at 
both  low  and  high  latitudes  and  inner  and  outer  heliospheric  locations  (Forsyth  1995:  Smith  et  al. 
1995;  McDonald  &  Lai  1995),  an  optimized  working  set  of  solar-minimum  solar-wind  turbulence 
and  HMF  parameters  can  be  collected. 

4.2.2  Modification  to  the  HMF  in  the  Polar  Regions 

This  modification  was  first  suggested  by  .lokipii  &  Kota  (1989)  to  allow  the  HMF  to  become 
larger  than  the  standard  Parker's  HMF  model  at  high  latitudes.  It  wa.s  motivated  by  recent 
observations  of  the  mean  and  fluctuating  components  of  the  HMF  at  high  latit  iirh's  (e.g..  Smith 
Bieber.  1991). 

The  modification  simply  amounts  to  adding  a  o-term  of  the  HMF  that  is  oc  r  and  implies 
a  non-zero  azimuthal  component  for  the  HMF  as  ^  >  0  (note  that  the  standard  Parker’s  spiral 
expression  suggests  that  B4,  0  as  0  0).  Smith  &  Bieber  (1991)  argued  that  one  could  invoke 

a  differential  solar  rotation  as  function  of  0  to  justify  this  modification. 

In  principle,  such  a  modification  is  straightforward  to  implement  in  a  numerical  code.  In  this 
code,  however,  and  since  we  calculate  all  diffusion  and  drift  terms  (which  are  affected  by  this 
modification)  analytically  so  as  to  minimize  the  numerical  error,  implementing  the  modification  is 
still  straightforward  but  somewhat  labor  (algebra)  intensive.  Finally,  we  note  that  this  particular 
modification  does  keep  the  large-scale  HMF  a  divergence  free  field,  with  or  w'ithout  the  presence  of 
the  sheet,  as  well  as  keep  V  •  (1^)  =  0.'  On  the  other  hand,  it  is  far  from  clear  that  this  particular 
modification  is  of  an.v  more  consequence  than,  .say.  the  (still  poorly  understood)  role  of  the  cross- 
field  diffusion  in  shaping  the  overall  transport  picture  of  galactic  cosmic  rays  in  a  three-dimensional 
heliosphere. 

4.2.3  The  Anomalous  CR  Component 

The  extension  of  the  code  to  include  the  anomalous  component  is  also,  in  principle,  straightfor¬ 
ward  since  the  transport  equation  for  this  component  is  essentially  identical  to  Eq.  (1)  except  for 
source  and  sink  terms  to  reflect  gains  and  losses  caiise<l  by  ionization,  hi  addition  to  these  terms, 
one  has  to  assume  an  initial  (or  injection)  distribution  function  for  the  .^CR  particle  in  the  vicinity 
of  the  shock  region  —around  80  .4U—  (Jokipii  1990). 


'Since  tliis  new  azimuthal  term  is  only  a  function  of  r.  the  overall  divergence  of  the  HMF  is  .still  zero.  ,Mso.  owing 
ic)  the  identity  V  •  V  X  [vector]  =  0.  and  from  Eq.  (T)  it  follows  that  (V'j)  is  still  divergence  li-ee  as  well. 
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Table  1.  Salient  physical  parameters  with  their  values  as  used  in  the  sample  calculations. 


PARAMETER 

VALUE 

UNITS 

TRANSPORTED  PARTICLE  Z-NUMBER 

1 

- 

TRANSPORTED  PARTICLE  A-NUMBER 

1 

— 

MAX.  PARTICLE  KINETIC  ENERGY 

1.000  X  10+- 

GeV/Nucl. 

MIN.  PARTICLE  KINETIC  ENERGY 

]. 000.x  lO"- 

GeV/NucI. 

INNER  HELIOSPHERIC  BOUNDARY 

1.000  X  10-- 

AU 

OUTER  HELIOSPHERIC  BOUNDARY 

1.000  X  10+* 

AU 

LISM-SPECTRUM  EXPONENT 

-2.65  X  lO^’ 

— 

STRENGTH  OF  DIFFUSION  TENSOR 

2.500  X  10+3 

AU'/yr 

PERP.  to  PARALLEL  DIFFUSION 

1.000  X  10-^ 

- 

STRENGTH  OF  DRIFT  VECTOR 

1.400  X  10+-' 

AU/yr 

STRENGTH  OF  HMF 

3.322  X  10+^ 

//G-AF- 

TILT  ANGLE  OF  SHEET  RMIN 

20 

des 

d'B-SPECTRUM  EXPONENT 

. .1.00  X  ID® 

- 

SOLAR-WIND  SPEED  'U  O-LAT. 

9.850  X  10+' 

AU/yr 

SOLAR-WIND  HELIOMAG-LAT  PAR. 

6.S00  X  10“^ 

- 

SOLAR  ROTATIONAL  SPEED 

9.450  X  10+* 

rad/yr 

BOUNDARY  CONDITION  (§■  RMIN 

F  =  0 

1/yr 

GRID  SIZE  [?•.</.  0,T] 

15  X  15  X  15  X  15 

- 
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A  perspective  of  the  magnetic  ciinent  sheet,  with  O'  given  by  Eip  (6).  in  cai  icsinn  coordinates. 
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Fig.  1.  Peclet  spectrum  for  H  as  a  function  of  kinetic  energy  and  polar 
angle  at  1  AU,  for  the  $  groups  ccdculated  according  to  (20). 


Fig.  2.  Peclet  spectrum  for  H  as  a  function  of  kinetic  energy  and  polar 
tingle  at  1  AU,  for  the  r  groups  calculated  according  to  (21). 
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10^  10^  10^  10^  10^ 
MeV/Nucleon 


Fig.  3.  Model  caloilations  of  solor-inininium  O-averaged  GCR-H  spectrum  at  ^  =  90 
deg  cuid  1  AU  (solid  curve)  and  at  100  Al’  (diisiied  curve).  Open  triangles  are  1987  solar- 
minimum  data  while  open  squares  are  1977  data.  Not  e  though  that  the  ccilcuiation.s  were 
performed  for  qA  >  0  .  i.e.,  1977  solar-minimum. 
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Fig.  4.  Calculated  ^averaged  solar-minimum  H  flu.x  at  1  AU  as  a  function  of  poltir 
angle  for  three  different  energies. 
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Local  ^  —  Gradient 


MeV/Nucleon 

Fig.  5.  Calculated  local  ^-gradient  as  a  function  of  energy.  ln[j(5‘’)/j(90‘^)]/85'^.  in 
%/deg.  at  1  AU. 
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Locail  1 — Gradient 


MeV/Nucleon 

FiR.  C.  Calculated  local  /-gradient.  lll[j(/’,v)/j(fx)]/(’’x.  ~  ''o)-  ;>•  runciioii  of 

energy,  in  %/AU,  at  i)  —  90  deg. 
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Appendix  1.  List  of  Main  FORTRAN??  Routines  for  code  PARADIGM 


PROGRAM  PARADIGM 

C .  .  . 

c . ; . . 

c. .  . 

C...  DYNAMIC  SOLAR -MODULATION  OF  GCR  IN  POLAR  SPHERICAL  COORDINATES 
C...  MODELLED  BY  A.  F.  BARGHOUTY 

C...  PHYSICS  DEPT.,  ROANOKE  COLLEGE,  SALEM,  VA  24153 

C...  [FOR  THE  COSMIC-RAY  SECTION,  NRL,  WASHINGTON,  DC] 

C...  [UNDER  NRL  GRANT  No.  N00014-95-1-G037] 

C.  .  . 

C. . .  June  1997 

C.  .  . 

C . 

c .  .  . 

C.  .  .  FLOWCHART  OF  THE  COMPLETE  SYSTEM  ORGANIZATION  AND  OPERATION  FOR  THE 
C...  SOLUTION  OF  THE  TRANSPORT  EQUATION  USING  THE  NUMERICAL  METHOD  OF 
C...  LINES  [cf.  Schiesser,  W.  E.  (1991)] 

U  .  .  . 


.BEGIN  EXECUTION  OF  MAIN. 
PROGRAM  PARADIGM  . 


C.  .  . 
C.  .  . 


C.  .  . 

C.  .  . 

C...  +  {+  DENOTES  AN  INPUT) 

C.  .  .  . 

C .  .  .  .  READ  THREE  . 

C.  .  .  . +  .DATA  LINES  (1)  . 

C. . .  .  . 

c. .  . 
c. .  . 
c. .  . 


c. 


.END  OF  RUNS . +  .  STOP  . 

.  LINE  READ  .  YES  . 


NO 


.PRINT  DATA. 
SUMMARY  . 


.  INCREMENT  THE  .  .  CALL  SUBROUTINE  INITIAL  . 

.RUN  COUNTER  (5).  .TO  INITIALIZE  THE  MODEL  (2). 


•PRINT  ERROR  SUMMARY. 
.  IF  REQUESTED  (4)  •  . 


YES  .  END  OF  RUN 

+ . (FINAL  VALUE 

.OF  TIME)  (3) 


NO 


COMMON/Y/  .  . 

+ .  CALL  SUBROUTINE  LSODES 

.  TO  INTEGRATE  THE  MODEL 

.  DIFFERENTIAL  EQUATIONS 

. +.OVER  ONE  PRINT  INTERVAL 

COMMON/ F/  . 


CALL  SUBROUTINE 
PDE  TO  COMPUTE 
THE  MODEL  TEMPORAL 
DERIVATIVES  (6) 


+ 


CALL  SUBROUTINE  PLOTS 
VIA  SUBROUTINE  PRINT 
TO  PRINT  THE  ENTIRE 
SOLUTION  VS  TIME 


[OPTIONAL] 


+ 


CALL  SUBROUTINE  PRINT 
TO  PRINT  THE  NUMERICAL 
SOLUTION  (INITIAL  CONDITIONS 
FOR  FIRST  CALL) .  STORE 
SOLUTION  FOR  SUBSEQUENT 
PLOTTING 


CALL  SYSTEM  UTILITIES 
TO  ASSIST  IN 
COMPUTATION  OF 
TEMPORAL  DERIVATIVES 
(7) 


EXPLANATORY  NOTES  FOR  THE  ABOVE  FLOWCHART 
(1)  LINE  1  -  TITLE(20)  (READ  VIA  900  FORMAT ( 2 0A4 ) ) 
LINE  2  -  T0,TF,TP  (READ  VIA  901  FORMAT ( 3E10 . 0 ) ) 


p  p  o  JO  f  p  o  /]  o  ^  p  r  p  o  jof  pn  Joi  pr  .  Jin  Jor  ho  !o'  ho 


C 


LINE  3  -  N,NMAX,NTYPE,NPRINT, IRRTYP, ERROR 

(READ  VIA  902  FORMAT (4 15 , 2X, 3A1 , ElO . 0 ) ) 


IF  **END  OF  RUNS**  IS  ENTERED  IN  COLUMNS  1  TO  11  OF  LINE  1  IN 
ANY  SET  OF  THREE  DATA  LINES,  PROGRAM  EXECUTION  IS  TERMINATED 
AND  LINES  2  AND  3  OF  THAT  SET  ARE  NOT  REQUIRED.  MULTIPLE 
SETS  OF  DATA  LINES  MAY  BE  USED,  THREE  LINES  PER  SET.  THE 
MAIN  PROGRAM  WILL  READ  EACH  SET  AND  EXECUTE  A  RUN  UNTIL  AN 
**END  OF  RUNS**  LINE  IS  READ. 


(2)  SUBROUTINE  INITIAL  IS  CALLED  ONCE  PER  RUN.  THEREFORE  DATA 
LINES  MAY  BE  READ  FROM  THIS  SUBROUTINE  TO  DEFINE  INITIAL 
PARAMETERS  OF  THE  MODEL  EQUATIONS  FOR  EACH  RUN.  THE  ADDI¬ 
TIONAL  DATA  LINES  WOULD  BE  PLACED  BEHIND  THE  THREE  BASIC 
DATA  LINES  OF  (1)  ABOVE. 

(3)  THE  END  OF  RUN  CONDITION  IS  T  GE  TF  WHERE  T  IS  THE  FIRST  ELE¬ 
MENT  IN  COMMON/T/  (GENERATED  BY  MAIN  PROGRAM  PARADIGM)  AND  TF 
IS  READ  FROM  DATA  LINE  (2)  OF  (1)  ABOVE. 

(4)  NPRINT  =  1  WILL  PRINT  A  SUMMARY  OF  THE  DEPENDENT  VARIABLES  IN 
COMMON/Y/  FOR  WHICH  THE  ESTIMATED  TEMPORAL  INTEGRATION  (TRUN¬ 
CATION)  ERROR  EXCEEDED  THE  MAXIMUM  PERMISSIBLE  VALUE,  ERROR, 
(READ  FROM  DATA  LINE  (3)  OF  (1)  ABOVE)  AT  ANY  POINT  DURING 
THE  RUN.  IF  NPRINT  =  0,  TEMPORAL  INTEGRATION  ERRORS  WILL  NOT 
BE  REPORTED. 

(5)  THE  RUN  COUNTER,  SET  BY  MAIN  PROGRAM  PARADIGM,  IS  THE  THIRD 
ELEMENT  IN  COMMON/T/  E.G.,  COMMON/ T / T, NF IN, NORUN 

(6)  THE  FUNDAMENTAL  LINKAGE  IN  THIS  SYSTEM  IS  THROUGH  COMMON/Y/ 
WHICH  CONTAINS  THE  MODEL  DEPENDENT  VARIABLE  VECTOR  AND  COMMON 
/F/  WHICH  CONTAINS  THE  VECTOR  OF  TEMPORAL  DERIVATIVES  OF  THE 
DEPENDENT  VARIABLE  VECTOR.  FOR  EXAMPLE,  THIS  LINKAGE  COULD 
BE  PROGRAMMED  AS 

COMMON/T/T, NF IN, NORUN/ Y/U(NEQ) /F/PUPT (NEQ) 

WHERE  THE  DEPENDENT  VARIABLE  VECTOR  U(NEQ)  IS  GENERATED  BY  THE 
TEMPORAL  INTEGRATOR,  SUBROUTINE  LSODES,  FROM  THE  DERIVATIVE 
VECTOR,  PUPT(NEQ),  GENERATED  BY  SUBROUTINE  PDE . 

MAIN  PROGRAM  PARADIGM  IS  THE  CALLING  PROGRAM  FOR  A  SERIES  OF  SUB¬ 
ROUTINES  WHICH  DEFINE  AND  INTEGRATE  THE  TEMPORAL  DIFFERENTIAL 

EQUATIONS.  THE  COMPLETE  PROGRAM  CONSISTS  OF  THE  FOLLOWING 

COMPONENTS 

(1)  MAIN  PROGRAM  PARADIGM  -  PERFORMS  OVERALL  CONTROL  OF  THE 
THE  TOTAL  PROGRAM. 


(2)  SUBROUTINE  INITIAL  -  SETS  THE  INITIAL  CONDITIONS  FOR  THE 
TEMPORAL  INTEGRATION. 

(3)  SUBROUTINE  PDE  -  DEFINES  THE  TEMPORAL  DERIVATIVE  VECTOR 
(PROVIDED  BY  THE  USER) .  ALSO  USES  SUPPLIED  FUNCTIONS. 


(4)  SUBROUTINE  DIFF  -  PERFORMS  SPATIAL  DIFFERENTIATION. 


C 


(5)  SUBROUTINE  LSODES  -  PERFORMS  THE  CENTRALIZED  TEMPORAL 
INTEGRATION. 


UU  UU  UU  UU  UU  UUhUU  UU  UU  CNUUO^UUmUU  uu 


c. .  . 


c. .  . 
c. .  . 
c. .  . 


COMMON/T/ 

T , TO , TF , NSTOP , NORUN 

/Y/ 

Y(IOOOO) 

/F/ 

F(IOOOO) 

COMMON/ 10/ 

NI,  NO 

DIMENSION  YV(IOOOO) ,  RWORK (4500000 ),  IWORK (10000 ) 

CHARACTER* 9  REALTIME , REALDATE 
EXTERNAL  FCN,  JAC 


ARRAY  FOR  THE  TITLE  (FIRST  LINE  OF  DATA) ,  CHARACTERS  END  OF  RUNS 
CHARACTER  TITLE (20) *4,  ENDRUN(3)*4 

VARIABLE  FOR  THE  TYPE  OF  ERROR  CRITERION 
CHARACTER* 3  ABSREL 


DEFINE  THE  CHARACTERS  END  OF  RUNS 
DATA  ENDRUN/'END  ','OF  RU','NS  '/ 

DEFINE  THE  INPUT/OUTPUT  UNIT  NUMBERS  AND  FILES 

NI=5 

NO=6 

CALL  TIME (REALTIME) 

CALL  DATE (REALDATE) 

OPEN(NI, FILE=' INPUT.DAT' , STATUS='OLD' ) 

OPEN (NO, FILE=' OUTPUT. DAT' , STATUS= ' NEW' ) 

INITIALIZE  THE  RUN  COUNTER 
NORUN=0 

BEGIN  A  RUN 
NORUN=NORUN+l 

INITIALIZE  THE  RUN  TERMINATION  VARIABLE 
NSTOP=0 

READ  THE  FIRST  LINE  OF  DATA 

READ (NI, 1000, END=999) (TITLE (I) , 1=1,20) 

TEST  FOR  END  OF  RUNS  IN  THE  DATA 
DO  2  1=1,3 

IF(TITLE(I) .NE.ENDRUN(I) )GO  TO  3 
CONTINUE 


AN  END  OF  RUNS  HAS  BEEN  READ,  SO  TERMINATE  EXECUTION 
STOP 

READ  THE  SECOND  LINE  OF  DATA 
READ(NI, 1001,END=999)T0,TF,TP 

READ  THE  THIRD  LINE  OF  DATA 

READ (NI, 1002,END=999)NEQN,NMAX, INT, I ERROR, ABSREL, ERROR 
PRINT  A  DATA  SUMMARY 

WRITE(NO, 1003)REALTIME, REALDATE, NORUN, (TITLE(I) ,1=1,20) , 

1  T0,TF,TP, 

2  NEQN , NMAX , INT , lERROR , ABSREL , ERROR 


INITIALIZE  TIME 
T=TO 


C 


C 


3. 


SET  THE  INITIAL  CONDITIONS 
CALL  INITIAL 

PRINT  THE  INITIAL  CONDITIONS 
CALL  PRINT (NI, NO) 

SET  THE  INITIAL  CONDITIONS  FOR  SUBROUTINE  LSODES 
DO  5  I=1,NEQN 
YV(I)  =Y(I) 

CONTINUE 

SET  THE  PARAMETERS  FOR  SUBROUTINE  LSODES 

TV=T0 

ITOL=l 

RTOL=ERROR 

ATOL=ERROR 

LRW=4500000 

LIW=10000 

IOPT=0 

ITASK=1 

I STATE =1 

MF=222 


INITIATE  THE  INTEGRATION 
TOUT=TV+TP 

CALL  SUBROUTINE  LSODES  TO  COVER  ONE  PRINT  INTERVAL 
CALL  LSODES { FCN , NEQN , YV , TV , TOUT , I TOL , RTOL , ATOL , ITASK , I STATE , 
1  I OPT , RWORK , LRW , I WORK , LIW , JAC , MF ) 

PRINT  THE  SOLUTION 
T=TV 

DO  6  1=1, NEQN 
Y { I ) =YV ( I ) 

CONTINUE 

CALL  PRINT (NI, NO) 


'i 

CL 


TEST  FOR  AN  ERROR  CONDITION 
IF ( ISTATE . LT . 0 ) THEN 

PRINT  A  MESSAGE  INDICATING  AN  ERROR  CONDITION 
WRITE (NO, 1004) ISTATE 

GO  ON  TO  THE  NEXT  RUN 
GO  TO  1 
END  IF 


'  . . .  CHECK  FOR  A  RUN  TERMINATION 
IF(NSTOP.NE.O)GO  TO  1 

.  . 

y  . .  THE  INTEGRATION  HAS  PROCEEDED  SATISFACTORILY,  SO  PREPARE  FOR  THE 
C. . .  NEXT  INTERVAL  IN  T .  NOTE  THAT  THE  FOLLOWING  CALCULATION  OF  THE 
C..  .  .  OUTPUT  TIME,  TOUT,  PRODUCES  THREE  OUTPUT  POINTS  FOR  EACH  DECADE 
...  IN  T 

TOUT=TV* (l.OE+01) ** (1 . OE+00/3 .E+00) 

C.  .  . 


C . . .  CHECK  FOR  THE  END  OF  THE  RUN 

IF(TV.LT. (TF-0.5E+00*TP) )GO  TO  4 

C.  .  . 

C...  THE  CURRENT  RUN  IS  COMPLETE,  SO  PRINT  THE  COMPUTATIONAL  STAT- 
C. . .  I ST ICS  FOR  LSODES  AND  GO  ON  TO  THE  NEXT  RUN 
CALL  TIME (REALTIME) 

CALL  DATE (REALDATE) 

WRITE (NO, 8)RWORK(ll) , IWORK(14) , IWORK(ll) , IWORK(12) , IWORK(13) , 

1  IWORK (17), IWORK (18), METHOD 

8  FORMAT (//I Ox, 'COMPUTATIONAL  STATISTICS  -  TIME  INTEGRATION:',/ 

2  lOxi'LAST  STEP  SIZE  =  ',1PE10.4,/ 

3  10x,'LAST  ORDER  OF  THE  METHOD  =  ',110,/ 


4 

lOx, 

'TOTAL  NUMBER  OF  STEPS  TAKEN 

= 

' , 110, / 

5 

lOx, 

'NUMBER 

OF 

FUNCTION  EVALUATIONS 

= 

MIO,/ 

6 

lOx, 

'NUMBER 

OF 

JACOBIAN  EVALUATIONS 

= 

' , 110, / 

7 

lOX, 

' LENGTH 

OF 

ARRAY [RWORK]  REQUIRED 

= 

' , 110, / 

8 

lOX, 

' LENGTH 

OF 

ARRAY [IWORK]  REQUIRED 

= 

' , 110, / 

9 

lOX, 

' INTEGRATION  METHOD 

' , 110, / 

1 

lOX, 

/ 

=  =  =  : 

- - - - - 

:  =  = 

81 

WRITE (NO , 81 ) NORUN, REALTIME , REALDATE 
FORMAT (1 OX, 'END  OF  RUN  NO.  -  ',12,' 

@  ',A8,2X,A9/) 

C.  .  . 
C.  .  . 

GO  TO  1 

C.  .  . 
C.  .  . 

FORMATS 

L.  .  .  . 

1000 

FORMAT ( 2 0A4) 

1001 

FORMAT (3E10 .0) 

1002 

FORMAT(4I5,2X,A3,E10.0) 

1003 

FORMAT ( / 

1  2X, A8, 2X,  A9, /, 

1  IX,'  RUN  NO.  -  ' ,I2,2X,20A4,//, 

2  IX,'  INITIAL  TIME  -  ',lpE10.3,/, 

3  IX,'  FINAL  TIME  -  ',lpE10.3,/, 

4  IX,'  PRINT  TIME  -  ',lpE10.3,/, 

5  IX,'  NUMBER  OF  DIFFERENTIAL  EQUATIONS  -  ',15,/, 

6  IX,'  PRINT  INTERVAL/MINIMUM  INTEGRATION  INTERVAL  -  ',15,/, 

7  IX,'  INTEGRATION  ALGORITHM  -  ',12,'  -  LSODES  ',/, 

8  IX,'  INTEGRATION  ERROR  MESSAGES  -  ',11,/, 

9  IX, '  ERROR  CRITERION  -  ' , A3 , / , 

A  IX,'  MAXIMUM  INTEGRATION  ERROR  -  ',1PE10.3) 

1004  FORMAT ( IX ,/, '  ISTATE  =  ',13,/, 

1  '  INDICATING  AN  INTEGRATION  ERROR,  SO  THE  CURRENT  RUN'  ,/, 

2  '  IS  TERMINATED.  PLEASE  REFER  TO  THE  DOCUMENTATION  FOR'  ,/, 

3  '  SUBROUTINE' ,/, 2 5X, 'LSODES' ,/, 

4  '  FOR  AN  EXPLANATION  OF  THESE  ERROR  INDICATORS'  ) 

END 

SUBROUTINE  FCN ( NEQN , TV , YV , YDOT ) 


C.  .  . 

C...  SUBROUTINE  FCN  IS  AN  INTERFACE  ROUTINE  BETWEEN  SUBROUTINES  LSODES 
C . . .  AND  DERV 

C.  .  . 

C. . .  DSS/2  COMMON  AREA 

COMMON/T/  T, TO, TF,NSTOP, NORUN 

1  /Y/  Y(l) 

2  /F/  F(l) 

C.  .  . 

C. . .  VARIABLE  DIMENSION  THE  DEPENDENT  AND  DERIVATIVE  ARRAYS 


)o  r 


C.  .  . 


C.  .  . 

c. . . 
£..  .  . 


2 


C 


C 


C 


DIMENSION  YV(NEQN) ,  YDOT(NEQN) 

TRANSFER  THE  INDEPENDENT  VARIABLE,  DEPENDENT  VARIABLE  VECTOR 

FOR  USE  IN  SUBROUTINE  DERV 

T=TV 

DO  1  I=1,NEQN 
Y ( I ) =YV ( I ) 

CONTINUE 

EVALUATE  THE  DERIVATIVE  VECTOR 
CALL  PDE 

TRANSFER  THE  DERIVATIVE  VECTOR  FOR  USE  BY  SUBROUTINE  LSODES 

DO  2  I=1,NEQN 

YDOT(I)=F(I) 

CONTINUE 

RETURN 

END 

SUBROUTINE  JAC 

SUBROUTINE  JAC  IS  A  DUMMY  ROUTINE  TO  SATISFY  THE  LOADER  (SINCE  JAC 
IS  DECLARED  AS  AN  EXTERNAL  IN  THE  MAIN  PROGRAM  LSODES) .  JAC  IS  NOT 
ACTUALLY  CALLED  UNLESS  AN  OPTION  OF  LSODE  IS  SELECTED  WHICH  REQUIRES 
THE  USER  TO  PROVIDE  THE  ODE  ANALYTICAL  JACOBIAN  MATRIX.  THIS  IS 
USUALLY  NOT  THE  CASE  SINCE  FOR  MOST  PROBLEMS,  THE  JACOBIAN  MATRIX 
IS  CALCULATED  INTERNALLY  BY  LSODE  USING  FINITE  DIFFERENCES 
RETURN 
END 


CONSTANTS'  FILE 


C . 

c 

C  Constant 
C 

c . 

c . 

c 

C  TILTO 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C  BO 

C 

C 

C 

C 

C 

C 

c 

C  Ws 

C 

C 

C  Vsw_0 
C 
C 
C 

C  WO 

c 

c 

c 

c 

c 

c 

C  DO 

c 

c 

c 

C  ETA 

C 
C 
C 

C  VO 

c 

c 

C  EO 

C 

C 

c  z 

c 

c 

C  A 

c 

c 


Description 


Tilt  angle  of  the  magnetic  neutral  sheet  AT  the  Sun. 

TILTO  corresponds  to  a  minimium  tilt  of  10  deg  at 
Solar  minimum  and  a  maximum  tilt  of  30  deg  at 
Solar  maximum,  i.e.,  the  tilt  of  the  sheet  correlates 
with  the  11 -yr  sunspot ’ cycle . 

•kie'kirieltie'kir'kie  *  ★ 

*******  at  TIME=0  yr,  SOLAR  MINIMUM  IS  ASSUMED  ********* 
**★★*★★*★*★★*★★★*★*★★*******★★*★★★★★★★**★★'*'***★★★★*★*'★**■■*★ 

TILTO  is  taken  here  to  be  20.0  deg. 
in  radians; 

Strength  of  the  photospheric  magnetic  field  x  solar 
radius'*‘2  -  corresponds  to  a  field  strength  of 
50.  micro-Gauss  at  Earth's  orbit.  (Could  be  positive 
or  negative  depending  on  polarity  and  the  22 -yr  solar 
cycle!)  Assumed  positive  over  the  first  11-yr  half 
solar  cycle  and  negative  over  the  second, 
in  (micro-Gauss  x  AU^2) ; 

Solar  rotational  speed  (corresponds  to  3.0E-6  rad/s) 
in  rad/yr; 

Solar-wind  speed  (assumed  radial  of  400  Km/s)  at  0-deg 
heliomagnet ic  latitude;  Vsw  has  a  heliomagnet ic  profile 
in  AU/yr; 

The  factor  WO  is  a  phenomenological  solar-activity 
factor  -  when  W0=0  the  solar  wind  is  latitude-indep . ; 
it  is  about  0.68  under  solar-minumum  conditions  for 
a  velocity  profile  that  increases  with  heliomagnetic 
latitude 
dimensionless ; 

Strength  of  local  diffusion  terms 
(corresponds  to  1.5E+22  cm''2/s) 
in  AU''2/yr; 

Relative  strength  of  normal  diffusion  term  to  tangential 
term  (in  local  solar-wind  frame) .  ETA  is  typically  <<  1 
dimensionless ; 

Strength  of  drift -velocity  terms 
in  AU/yr 

amu  Rest -mass 
in  GeV/c''2; 

Atomic  number  of  particle  being  transported 
dimensionless; 

Mass  number  of  particle  being  transported 
dimensionless ; 


uq  'U|  uLl  jul  oU  juI  u(  ud  jul  ud  jul  d-  uL  ud  'bul  ju 


Rtnax 


Upper  range  of  radial  distance  from  the  SUN,  i.e., 

C  outer  Heliospheric  boundary 

in  AUs 

Rmin  Lower  range  of  radial  distance  from  the  SUN,  i.e., 

inner  Heliospheric  boundary  -  nonzero 
in  AUs 

Pmax  Upper  range  of  polar  angle  -w.r.t.  rotational  axis  of 

the  SUN 
in  rad 

Pmin  Lower  range  of  polar  angle  -w.r.t.  rotational  axis  of 

the  SUN 
in  rad 

Amax  Upper  range  of  azimuthal  angle  around  the  Sun 

in  rad 

Amin  Lower  range  of  azimuthal  angle  around  the  Sun 

in  rad 

Emax  Upper  range  of  particle's  kinetic  energy 

in  GeV/nucleon 

Emin  Lower  range  of  particle's  kinetic  energy 

in  GeV/nucleon 

EXPONENT  This  is  the  exponent  in  total  energy  of  the  assumed 
LISM  density  -  typically  -2.65  to  -2.75 

This  is  the  exponent  in  the  power  spectrum  of  the  fluctuating 
component  of  the  Heliospheric  magnetic  field  -  a  purely 
Kolmogorov  spectrum  has  EXP=5/3;  a  Kraichnan  spectrum  has 
EXP=3/2,  etc. 

This  is  the  assumed  boundary  condition  at  R=RMIN;  when 
BC='U0'  the  model  assumes  U(RMIN)=0.  for  all  times,  i.e., 
the  density  function  and  hence  the  intensity  is  assumed 
zero  at  R=RMIN,  and  when  BC='SF',  the  model  assumes  the 
streaming  flux  being  zero  at  R=RMIN,  i.e.,  absorbing  Sun, 
which  satisfies  Liouville's  theorem. 

Speed  of  light  in  AU/yr 


Grid:  3 -dimensional  heliocentric  polar  spherical 

coordinate  system)  +  Energy 
(  Radial,  POLAR,  Azimuthal,  Kinetic  Energy) 
PARAMETER (  NR=15,  NP=15,  NA=15,  NE=15  ) 

CHARACTER* 2  BC 

LISM-spectrum: 

DATA  EXPONENT/ -2 . 65/ 

Fluctuating  component  of  HMF : 

DATA  EXP/1.000/ 

Regular  component  of  HMF : 

DATA  B0,TILT0/33 .218, .349/ 

Solar  Wind: 

DATA  WS,Vsw  0,W0/94.5,84.4,.68/ 


Drift : 


no  no  on  n  on  on  on 


DATA  V0/1.4E+4/ 

Diffusion : 

DATA  DO , ETA/2 . 5E+3 , .10/ 

Particle : 

DATA  E0,ZNo,ANo/0.938,l. ,1./ 

3-D  Grid  Ranges: 

DATA  RMAX, RMIN/100 . , . 01/ 

DATA  PMAX,PMIN/l.57080,8.72664E-2/ 
DATA  AMAX, AMIN/3 .14159, .0/ 

Energy  Range : 

DATA  EMAX, EMIN/100 . , .01/ 

Boundary  condition  at  R=RMIN: 

DATA  BC/'SF'/ 

Physical  and  mathematical  constants 
DATA  c, PI/6. 31E+4, 3. 14159/ 


'EST  RUN(S) 
O.E+00 
-^06251000 
ND  OF  RUNS 


FOR  THE  SOLUTION  OF  THE  JP  EQUATION 
.OOE+00  l.E-6 

2  1  REL  l.E-02 


c 

c 


COMMON  BLOCKS  FILE 


COMMON/ T/ 
/Y/ 
/F/ 
/S/ 


T,  TO,  TF,  NSTOP,  NORUN 
U{NR,NP,NA,NE) , TYME ( 0 : 1000000 ) 
UT{NR,NP,NA,NE) 


UR(NR,NP,NA,NE) 
UP (NR,NP,NA,NE) 
UA(NR,NP,NA,NE) 
AD(NR,NP,NA,NE) 
UO (NR,NP,NA,NE} 


URA(NR,NP,NA,NE) , 


ADER(NR,NP,NA,NE) , 


URR{NR,NP,NA,NE) 

UPP {NR,NP,NA,NE) 

UAA(NR,NP,NA,NE) 

ADE (NR,NP,NA,NE) 

TLUX_r (NR , NP , NA, NE ) , 

FLUX_p (NR,NP,NA,NE) ,  FLUX_a (NR, NP, NA, NE) , 

ANISOT_r (NR,NP,NA,NE) ,  ANISOT_p (NR,NP,NA,NE) , 
ANISOT_a(NR,NP,NA,NE) ,  FLUX (NR, NP, NA, NE) 

/C/  RMIN,  RMAX,  PMIN,  PMAX,  AMIN,  AMAX,  EMIN,  EMAX, 
R(NR),  P(NP),  A{NA),  E(NE),  lEO,  JEO,  KEO,  BC 
COMMON/PDE/DIFFUSION, CONVECTION, DRIFT, ADIABATIC, SOURCE , 

1  R_DIFFUSION, A_DIFFUSION, P_DIFFUSION, R_DRIFT, P_DRIFT, 

2  A  DRIFT 


SUBROUTINE  INITIAL 


INCLUDE  'CONSTS. MODEL' 

INCLUDE  'COMMON. MODEL' 

c 

__  common  nprint ,  ip 

COMMON/ 10/  NI,NO 

C.  .  . 

—  NPRINT=1 

C. . .  COMPUTE  THE  RADIAL,  POLAR,  AZIMUTHAL,  AND  ENERGY  POSITIONS: 

.  NOTE:  RADIAL  AND  ENERGY  POSITIONS  ARE  LOGARITHMICALLY  SPACED! 

DELTA_R=(RMAX/RMIN)**<1./(NR-1.) ) 

IE0=1 
DO  1  =  1,  NR 

IF(I.EQ.l)  THEN 
R(l) =RMIN 
~  ELSE 

R ( I ) =R ( I - 1 ) *DELTA_R 
END  IF 

_  IF  (R(I) .GT. 0.75. AND. R(I) .LT. 1.15)  IE0=I 

END  DO 

JEO=l 

—  DO  J=1,NP 

P  ( J)  =  ( J-1)  *  (PMAX-PMIN)  /  (NP-1)  +PMIN 
END  DO 

£..  .  . 

KEO=l 
DO  K=1,NA 

A(K) = (K-1) * (AMAX-AMIN) / (NA-1) +AMIN 
END  DO 

DELTA_E= (EMAX/EMIN) **(!./ (NE-1.  )  ) 

—  DO  L=1,NE 

IF(L.EQ.l)  THEN 
E (1) =EMIN 
_  ELSE 

E (L)  =E (L-1) *DELTA_E 
END  IF 
END  DO 

\..  SET  THE  INITIAL  CONDITIONS  AT  TIME=0,  I.E.,  LOCAL- ISM  SPECTRUM: 

DO  1=1, NR 
_  DO  J=1,NP 

DO  K=1,NA 
DO  L=1,NE 

_  IF(I.EQ.NR)  THEN 

UO (I, J,K,L) =SPECTRUM(E (L) ) *VELOCITY(E (L) ) / (4 . *PI) 

U(I,  J,K,L)=U0  (I,  J,K,L) 

ELSE 

—  U(I,  J,K,L)  =0. 

END  IF 

^  END  DO 

_  END  DO 

END  DO 
END  DO 

C.  .  . 


RETURN 
END  . 


SUBROUTINE  PDE 


LOGICAL  EARTH_ORBIT 
INCLUDE  ' CONSTS. MODEL' 

INCLUDE  'COMMON. MODEL' 

ic=ic+l  ! COUNTER 
TYME (IC) =T 

ZERO  ALL  TERMS; 

U  is  omni-directional,  differential  CR-intensity  function  that  is 
being  transported; 

[in  no.  of  particles/ (m'^2-sr-s-GeV/N)  ]  : 

DO  1=1, NR 
DO  J=1,NP 
DO  K=1,NA 
DO  L=1,NE 

UT(I,  J,K,L)  =0  . 

UR(I, J,K,L) =0. 

URR(I, J,K,L) =0. 

URA(I, J,K,L) =0. 

UP(I,  J,K,L)=0. 

UPP(I,  J,K,L)  =0. 

UA(I,  J,K,L)  =0. 

UAA(I, J,K,L) =0. 

AD(I,  J,K,L)=0. 

ADE(I,  J,K,L)  =0. 

ADERd,  J,K,L)  =0  . 

R_DIFFUSION=0 . 

P_DIFFUSION=0. 

A_DIFFUSION=0 . 

DIFFUSION=0 . 

CONVECTION=0 . 

R_DRIFT=0 . 

P_DRIFT=0 . 

A_DRIFT=0. 

DRIFT=0 . 

ADIABATIC=0 . 

END  DO 
END  DO 
END  DO 
END  DO 

SET  BOUNDARY  CONDITION  AT  R=RMAX : 

DO  J=1,NP 
DO  K=1,NA 
DO  L=1,NE 

U(NR, J,K,L) =U0 (NR, J,K,L) 

END  DO 
END  DO 
END  DO 

SET  THE  BOUNDARY  CONDITION  AT  E=EMAX 
DO  1=1, NR 
DO  J=1,NP 
DO  K=1,NA 

U(I, J,K,NE) =U0 (NR, J,K,NE) 

END  DO 
END  DO 
END  DO 


no  noooooonn  no  on 


C. . .  SET -BOUNDARY  CONDITION  AT  R=RMIN  IF  BC='U0' : 

IF{BC.EQ. 'UO' )  THEN 

C. . .  NOTE;  This  BC  corresponds  to  U(RMIN)=0.  for  all  T! 

DO  J=1,NP 
DO  K=1,NA 
DO  L=1,NE 

U(l,  J,K,L)  =0. 

END  DO 
END  DO 
END  DO 
END  IF 

ADIABATIC- COOLING  TERM: 

DO  1=1, NR 
DO  J=1,NP 
DO  K=1,NA 
DO  L=1,NE 

AD(I,  J,K,L)  =ALPHA{E(L)  )  *E  (L)  *U  (I ,  J,  K,  L) 

END  DO 
END  DO 
END  DO 
END  DO 

COMPUTE  THE  FIRST-ORDER  DERIVATIVES  IN  E : 

CALL  DSS_4d (LOG (EMIN) ,LOG(EMAX) , NR, NP , NA, NE , 4 , AD, ADE , 0 . ) 

DO  1=1, NR 
DO  J=1,NP 
DO  K=1,NA 
DO  L=1,NE 

ADE(I,  J,K,L)=ADE(I,  J,K,L)  /E(L) 

END  DO 
END  DO 
END  DO 
END  DO 

COMPUTE  THE  SECOND-ORDER  MIXED  DERIVATIVES  IN  E-R  IF  BC='SF': 

IF (BC.EQ. ' SF' )  THEN 

CALL  DSS_4d(LOG(RMIN) ,LOG(RMAX) , NR, NP, NA, NE, 1 , ADE , ADER, 0 . ) 

This  mixed  derivative  is  needed  as  a  boundary  condition  at  R=RMIN 
when  the  "S"treaming  "F"lux  is  assumed  zero  at  R=RMIN,  rather  than 
density  U  being  assumed  zero.  "SF"  BC  satisfies  Liouville's  theorem, 
whereas  BC  "UO"  does  not! 

DO  1=1, NR 
DO  J=1,NP 
DO  K=1,NA 
DO  L=1,NE 

ADERd,  J,K,L)  =ADER(I,  J,K,L)  /R(I) 

END  DO 
END  DO 
END  DO 

END  DO  ^  - 

END  IF 

COMPUTE  THE  FIRST-ORDER  DERIVATIVES  IN  R  FOR  DIFFUSION: 

CALL  DSS_4d(LOG(RMIN) ,LOG(RMAX) ,NR,NP,NA,NE,1,U,UR,0.) 

DO  1=1, NR 


n  r 


DO  J=1,NP 
DO  K=1,NA 
DO  L=1,NE 

UR{I,  J,K,L)  =UR(I,  J,K,L)  /R(I) 

END  DO 
END  DO 
END  DO 
END  DO 

COMPUTE  THE  SECOND-ORDER  DERIVATIVES  IN  R  FOR  DIFFUSION: 

CALL  DSS_4d{LOG(RMIN) ,LOG{RMAX) , NR, NP, NA, NE , 1 , UR, URR, 0 . ) 

DO  1=1, NR 
DO  J=1,NP 
DO  K=1,NA 
DO  L=1,NE 

URR(I, J,K,L) =URR(I, J,K,L) /R(I) 

END  DO 
END  DO 
END  DO 
END  DO 

COMPUTE  THE  SECOND-ORDER  MIXED  DERIVATIVES  IN  R-A  FOR  DIFFUSION: 
CALL  DSS_4d(AMIN,AMAX,NR,NP,NA,NE,3,UR,URA,0. ) 

COMPUTE  THE  FIRST-ORDER  DERIVATIVES  IN  R  FOR  CONVECTION: 

CALL  DSS_4d(LOG(RMIN) ,LOG(RMAX) , NR, NP, NA, NE , 1 , U, UR, -1. ) 

DO  1=1, NR 
DO  J=1,NP 
DO  K=1,NA 
DO  L=1,NE 

UR(I,  J,K,L)  =UR(I,  J,K,L)  /R(I) 

END  DO 
END  DO 
END  DO 
END  DO 

COMPUTE  THE  FIRST-ORDER  DERIVATIVES  IN  P  FOR  DIFFUSION: 

CALL  DSS_4d(PMIN, PMAX, NR, NP , NA, NE , 2 , U, UP, 0 . ) 

SET  BOUNDARY  CONDITIONS  AT  P=PMIN  and  P=PMAX: 

DO  1=1, NR 
DO  K=1,NA 
DO  L=1,NE 

UP(I, 1,K,L) =0. 

UP(I,NP,K,L) =0. 

END  DO 
END  DO 
END  DO 

COMPUTE  THE  SECOND-ORDER  DERIVATIVES  IN  P  FOR  DIFFUSION: 

CALL  DSS_4d(PMIN,PMAX,NR,NP,NA,NE,2,UP,UPP, 0. ) 

COMPUTE  THE  FIRST-ORDER  DERIVATIVES  IN  P  FOR  CONVECTION: 

CALL  DSS_4d(PMIN,PMAX,NR,NP,NA,NE,2,U,UP, 0. ) 

SET  BOUNDARY  CONDITIONS  AT  P=PMIN  and  P=PMAX: 

DO  1=1, NR 
DO  K=1,NA 
DO  L=1,NE 

UP(I,1,K,L)=0. 


on  non  ononon  onon 


UP(I,NP,K,L)  =0. 

END  DO 
END  DO 
END  DO 

COMPUTE  THE  FIRST-ORDER  DERIVATIVES  IN  A  FOR  DIFFUSION: 
CALL  DSS_4d{AMIN,AMAX,NR,NP,NA,NE,3,U,UA,0. ) 

SET  BOUNDARY  CONDITIONS  AT  A=AMIN  and  A=AMAX: 

DO  1=1, NR 
DO  J=1,NP 
DO  L=1,NE 

UA(I, J, 1,L) =0. 

UA(I, J,NA,L) =0. 

END  DO 
END  DO 
END  DO 

COMPUTE  THE  SECOND-ORDER  DERIVATIVES  IN  A  FOR  DIFFUSION: 
CALL  DSS_4d(AMIN, AMAX,NR,NP,NA,NE, 3,UA,UAA, 0 . ) 

COMPUTE  THE  FIRST-ORDER  DERIVATIVES  IN  A  FOR  CONVECTION: 
CALL  DSS_4d(AMIN,AMAX,NR,NP,NA,NE,3,U,UA, 0 . ) 

SET  BOUNDARY  CONDITIONS  AT  A=AMIN  and  A=AMAX: 

DO  1=1, NR 
DO  J=1,NP 
DO  L=1,NE 

UA(I, J, 1,L) =0 . 

UA(I, J,NA,L) =0. 

END  DO 
END  DO 
END  DO 

ASSEMBLE  THE  PDE : 

DO  1=1, NR 
DO  J=1,NP 
DO  K=1,NA 
DO  L=1,NE 


C 

C.  .  . 
C.  .  . 
C.  .  . 

C 

C.  .  . 


The  collective  sum  of  radial  diffusion  terms: 


+ 

+ 

+ 

+ 


+ 


R_DIFFUSION=UR(I, J,K,L) * (2.*Drr (R(I) , P ( J) , A (K) , E (L) , T) /R ( I ) + 

Drr_r (R(I) , P ( J) , A (K) , E (L) , T) ) + 

URR(I, J,K,L) *Drr (R(I) , P ( J) , A (K) , E (L) , T) + 

(UA(I, J,K,L) * (Dra(R(I) , P ( J) , A (K) , E (L) , T) /R { I ) + 
Dra_r (R(I) ,P(J) , A (K) , E (L) , T) ) + 

URAd,  J,K,L)  *Dra(R{I)  ,P(J)  ,  A  (K)  ,  E  (L)  ,  T)  )  / 

(R(I)  *SIN(P(J)  )  ) 


NOTE:  URA ( I , J, K, L) =UAR (I , J, K,L) ;  the  term  appears  both 
in  the  radial  diffusion  terms  as  well  as  in  the 
azimuthal  diffusion  terms. 


The  collective  sum  of  polar  diffusion  terms: 
P_DIFFUSION=UP(I, J,K,L) / (R(I) **2) * 

+  (Dpp(R(I)  ,P(J)  ,A(K)  ,E(L)  ,T)  *COS(P(J)  )  / 

+  SIN(P(J)  ) +Dpp_p{R(I)  ,P(J)  ,A(K)  ,E(L)  ,T)  )  + 

+  UPPd,  J,K,L)  /  (R{I)  **2) 

+  *Dpp(Rd)  ,P(J)  ,A(K)  ,E(L)  ,T) 


|n  Tv,  p  n  1  o  c  h  o  V  >  o  lor  n  )  o  (  1  n  1  o  In  r 


+ 

+ 

+ 


c. .  . 


The  collective  sum  of  azimuthal  diffusion  terms: 

A_DIFFUSION=UAA(I, J,K,L) *Daa{R(I) ,P{J) ,A(K) ,E(L) ,T) / 
(R(I)*SIN(P(J)))**2+ 

URA(I, J,K,L) *Dar(R{I) , P ( J) , A (K) , E (L) , T) / 
(R(I)  *SIN(P(J)  )  ) 


The  collective  sum  of  all  diffusion  terms  in  3d: 

DIFFUSION=R  DIFFUSION+P  DIFFUSION+A  DIFFUSION 


The  convection  term  due  to  a  radial  solar-wind  velocity  profile: 
CONVECTION=UR(I, J,K,L) *Vsw(P(J) ) + 

+  2.*U(I,J,K,L)*Vsw(P(J)  )/R(I) 


C 


The  radial,  polar,  and  azimuthal  drift  terms: 

R_DRIFT=VDR(R(I) , P ( J) , A (K) , E (L) , T) *UR{I, J,K,L) 

P_DRIFT=VDP(R(I) ,P(J) ,A(K) ,E(L) ,T) *UP(I, J,K,L) /R(I) 

A_DRIFT=VDA(R(I) ,P(J) ,A(K) ,E{L) ,T) *UA(I, J,K,L) / 

+  (R(I) *SIN(P(J) ) ) 


The  collective  sum  of  drift  terms  in  3d: 
DRIFT=R  DRIFT+P  DRIFT+A  DRIFT 


The  adiabatic  cooling  term: 

ADIABATIC= (2 , *Vsw(P (J) ) / (3 . *R(I) ) ) *ADE (I, J,K,L) 


Assumed  sources  [apart  from  ISM  boundary  condition]  if  any: 
SOURCE=0 . 


Calculate  Streaming  Flux  Vector: 

.  ...  Radial  component  of  the  streaming  flux  vector: 

FLUX_r  (I,  J,K,L)  =-Drr  (R(I)  ,  P  ( J)  ,  A (K)  ,  E  (L)  ,  T)  *UR(I,  J,K,L)  - 
+  Dra(R(I)  ,  P  ( J)  ,  A  (K)  ,  E  (L)  ,  T)  *UP(I,  J,K,L)  /  (R(I)  *SIN(P(J)  )  )  + 

+  VDR(R(I)  ,P(J)  ,A(K)  ,E(L)  ,  T)  *U  ( I ,  J,  K,  L)  + 

+  Vsw(P(J)  )  *  (U(I,  J,K,L)  -  {l./3.)*ADE(I,  J,K,L)  ) 

.  ...  Polar  component  of  the  streaming  flux  vector: 

FLUX_p(I,  J,K,L)  =-Dpp(R(I)  ,  P  ( J)  ,  A  (K)  ,  E  (L)  ,  T)  *UP(I,  J,K,L)  /R(I)  + 

+  VDP(R(I)  ,P(J)  ,A(K)  ,E(L)  ,T)  *U(I,  J,K,L) 

.  ...  Azimuthal  component  of  the  streaming  flux  vector: 

FLUX_a(I,  J,K,L)=-Dar(R(I)  ,P(J)  ,A(K)  ,E{L)  ,  T)  *UR  ( I ,  J,  K,  L)  - 
+  Daa(R(I)  ,  P  ( J)  ,  A  (K)  ,  E  (L)  ,  T)  *UA ( I ,  J,  K,  L)  /  (R  ( I )  *SIN  (P  ( J)  )  )  + 

+  VDA(R(I)  ,P(J)  ,A{K)  ,E(L)  ,  T)  *U  ( I ,  J,  K,  L) 

Calculate  Anisotropy  Vector: 

.  . . .  Radial  component  of  the  anisotropy  vector: 

IF(U(I, J,K,L) .NE.O. )  THEN 

ANISOT_r (I, J,K,L) =3 . *FLUX_r { I , J, K, L) / (4 . *PI*U ( I ,  J,  K,  L)  ) 

.  . . .  Polar  component  of  the  anisotropy  vector: 

ANISOT_p(I, J,K,L) =3.*FLUX_p(I, J,K,L) / (4 . *PI*U ( I , J, K, L) ) 

.  . . .  Azimuthal  component  of  the  anisotropy  vector: 

ANISOT_a(I, J,K,L) =3 . *FLUX_a ( I , J, K, L) / (4 . *PI*U { I , J, K, L) ) 

END  IF 


C 


c 

C...  Assume  modulation  is  negligible  at  heliospheric  boundary, 
C...  and  at  the  highest  energy  [lOO  GeV/Nucl.]: 
IFd.EQ.l.AND.BC.EQ. 'UO' )  THEN 
U(I,J,K,L)=0. 

UT{I,  J,K,L)=0. 

ELSE  IFd.EQ.l.AND.BC.EQ. 'SF' )  THEN 

UTd,  J,K,L)  =-Vsw(P(J)  )  *ADERd,  J,K,L)  /3  . 

ELSE  IFd.EQ.NR)  THEN 

U  (I ,  J,  K,  L)  =U0  (NR,  J,  KT,  L) 

UTd,  J,K,L)=0. 

ELSE  IF(L.EQ.NE)  THEN 

Ud,  J,K,L)=U0  (NR,  J,K,L) 

UTd,  J,K,L)  =0. 

ELSE 

C. . .  JP  transport -equation  in  3d  at  this  T: 

C . 

UT (I , J, K, L) =DIFFUSION-CONVECTION-DRIFT+ADIABATIC+SOURCE 

C . 

END  IF 
C 

FLUXd,  J,K,L)  =Ud,  J,K,L) 

C 

c  CALL  WARNd,  J,K,L,  IC) 

c  CALL  DIAGd,  J,K,L,  IC) 

C 

END  DO 
END  DO 
END  DO 
END  DO 

C.  .  . 

RETURN 

END 


Gp=WO*SIN  (2  .  *POIiAR)  /  (1 .  +WO*SIN  (PI/2  .  -POLAR)  **2) 

Note  that  Gp=0  if  W0=0,  i.e.,  latitude- indep .  solar-velocity! 
Gpp=l . +RADIAL* (Ws/Vsw (POLAR) ) *TILT (TIME) * 

COS (AZIMUTH-RADIAL* (Ws/Vsw (POLAR) ) ) *Gp 

POLARO=Pl/2. +TILT (TIME) *SIN (AZIMUTH-RADIAL* (Ws/Vsw (POLAR) ) ) 

BO=SIGN(BO,COS(PI*TIME/22.) ) 

VD0=V0* (2 . *BETA*Pc) / (3 . *ZNo*B0* (l.+Gl**2) **2) *RADIAL 
VD0=VD0  *DIRAC ( POLAR , POLARO ) 

VDR=VDSR (RADIAL , POLAR , AZIMUTH , ENERGY , TIME ) 

*SHEET (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) +  VD0*Gl*Gpp 

RETURN 

END 

FUNCTION  VDSP (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

Polar  component  of  the  drift  velocity 
in  a  single  sector 

INCLUDE  'CONSTS. MODEL' 

GAMMA=1 . +ENERGY/E0 
BETA=SQRT ( 1 . -GAMMA* *  - 2 ) 

Pc=ANo*SQRT (ENERGY* (ENERGY+2 . *E0) ) 

G1=RADIAL* (Ws/Vsw (POLAR) ) *SIN (POLAR) 

G2=RADIAL* (Ws/Vsw (POLAR) )*COS( POLAR) 

B0=SIGN (BO , COS (PI*TIME/22 . ) ) 

VD0=V0* (2 . *BETA*Pc) / (3 . *ZNo*B0* (l.+Gl**2) **2) *RADIAL 

VDSP=VD0*G1* (2 .+G1**2) 

RETURN 

END 

FUNCTION  VDP (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

POLAR  component  of  the  drift  velocity 
given  the  neutral  magnetic  sheet  SHEET 

INCLUDE  ' CONSTS . MODEL ' 

GAMMA=1 . +ENERGY/E0 
BETA=SQRT ( 1 . -GAMMA* *  - 2 ) 

Pc=ANo* SORT (ENERGY* (ENERGY+2 . *E0) ) 

G1=RADIAL* (Ws/Vsw (POLAR) ) *SIN (POLAR) 

G2=RADIAL* (Ws/Vsw ( POLAR) ) *COS (POLAR) 

POLARO=Pl/2.+TILT(TIME) *SIN (AZIMUTH-RADIAL* (Ws/Vsw (POLAR) ) ) 
B0=SIGN (BO , COS (PI*TIME/22 . ) ) 

VD0=V0* (2 . *BETA*Pc) / (3 . *ZNo*BO* (l.+Gl**2) **2) *RADIAL 
VD0=VD0*DIRAC (POLAR, POLARO) 

VDP=VDSP (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

*SHEET (RADIAL, POLAR, AZIMUTH, ENERGY, TIME)  + 
VDO*RADIAL*TILT (TIME) *COS (AZIMUTH-RADIAL* (Ws/Vsw ( POLAR) ) ) 
Gl* (Ws/Vsw (POLAR) )* (Gl**-2-l. ) 


c 


c 

c. 

c. 

c 

c 


c 

c. 

c 


c 

c 


c 

c. 

c. 

c 

c 


c 

c. 


c 


c 


c 


c 

c. 

c 


RETURN 

END 

FUNCTION  VDS A  (RADIAL ,  POLAR ,  AZIMUTH ,  ENERGY ,  TIME )' 

Azimuthal  component  of  the  drift  velocity 
in  a  single  sector 

INCLUDE  ' CONSTS. MODEL' 

GAMMA=1 . +ENERGY/E0 
BETA=SQRT ( 1 . -GAMMA* * - 2 ) 

Pc=ANo*SQRT (ENERGY* (ENERGY+2 . *E0) ) 

G1=RADIAL* (Ws/Vsw( POLAR) ) *SIN (POLAR) 

G2=:RADIAL*  (Ws/Vsw  ( POLAR)  )*COS(  POLAR) 

G=W0* SIN (POLAR) **2/ (1 . +W0*SIN (PI/2 . -POLAR) **2) 

Note  that  G=0  if  W0=0,  i.e.,  latitude - indep .  solar- velocity ! 

BO=SIGN(BO, COS (PI*TIME/22 .  )  ) 

VD0=V0* (2 . *BETA*Pc) / (3 . *ZNo*B0* (l.+Gl**2) **2) *RADIAL 

VDSA=VD0*G1*G2* (1 . +2 . *G) 

RETURN 

END 

FUNCTION  VDA (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

Azimuthal  component  of  the  drift  velocity 
given  the  neutral  magnetic  sheet  SHEET 

INCLUDE  'CONSTS. MODEL' 

GAMMA=1 . +ENERGY/E0 
BETA=SQRT ( 1 . -GAMMA* * - 2 ) 

Pc=ANo*SQRT (ENERGY* (ENERGY+2 . *E0) ) 

G1=RADIAL* (Ws/Vsw (POLAR) ) *SIN (POLAR) 

G2=RADIAL* (Ws/Vsw ( POLAR) ) * COS (POLAR) 

Gp=WO*SIN(2.*POLAR) / ( 1 . +W0*SIN (PI/2 . -POLAR) **2) 

Note  that  Gp=0  if  W0=0,  i.e.,  latitude-indep .  solar-velocity 
Gpp=l . +RADIAL* (Ws/Vsw (POLAR) ) *  TILT (TIME) * 

+  COS (AZIMUTH -RADIAL* (Ws/Vsw (POLAR) ) ) *Gp 

POLARO=PI/2 . +TILT (TIME) *SIN (AZIMUTH -RADIAL* (Ws/Vsw ( POLAR) ) ) 

BO=SIGN(BO,COS(PI*TIME/22.) ) 

VD0=V0* (2 . *BETA*Pc) / (3 . *ZNo*BO* (1 . +G1**2 ) **2 ) *RADIAL 
VD0=VD0 *DIRAC ( POLAR , POLARO ) 

VDA=VDSA (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

+  *SHEET (RADIAL, POLAR,AZIMUTH, ENERGY, TIME)  +  VD0*Gpp 

RETURN 

END 

FUNCTION  VDS (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

Drift  speed  in  a  single  sector 

VDS =SQRT ( VDSR (RADIAL , POLAR , AZ IMUTH , ENERGY , T IME ) *  * 2 


FUNCTION  ALPHA (ENERGY) 


INCLUDE  'CONSTS. MODEL' 

ALPHA= (ENERGY+2 . *E0) / (ENERGY+EO) 

RETURN 

END 

FUNCTION  SPECTRUM (ENERGY) 

INCLUDE  ' CONSTS. MODEL' 

CHOSEN  SO  THAT  AT  100  GeV,  THE  LISM  UNMODULATED 
INTENSITY  MATCHES  THE  OBSERVED  INTENSITY,  WHERE  THE 
MODULATION  EFFECTS  ARE  ASSUMED  NEGLIGIBLE! 

CONST=1.23*2.*l.E-3/l.  ! f or  p;  lism-exp=-2 . 65 

CONST=l. 5852*1. 23*2. *l.E-3/l9. 018  ! f or  He 

CONST=l. 5852*1. 23*2. *l.E-3/3651.  !for  Fe 


SPECTRUM=CONST* (ENERGY+EO) **EXPONENT 

RETURN 

END 

FUNCT ION  Velocity( ENERGY ) 

To  convert  particle  density  to  particle  intensity: 

INCLUDE  'CONSTS. MODEL' 

GAMMA=1 . +ENERGY/E0 
BETA=SQRT (1 . -GAMMA** -2) 

Velocity=BETA*c 

RETURN 

END 

FUNCTION  Vsw( POLAR) 

Solar-wind  velocity  profile  -  assumed  radial  but  heliomagnet 
latitude  dependent: 

INCLUDE  ' CONSTS . MODEL ' 

Vsw=Vsw_0* (l.+W0*SIN( PI/2. -POLAR) **2) 

RETURN 

END 

FUNCTION  Br (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

Radial  component  of  the  solar  magnetic  field 
in  a  single  sector 

INCLUDE  ' CONSTS. MODEL' 

B0=SIGN(B0, COS (PI*TIME/22 . ) ) 

Br=B0 /RADIAL* *2 


RETURN 

END 

C 

FUNCTION  Bphi (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

C. . .  Azimuthal  component  of  the  solar  magnetic  field 
C. . .  in  a  single  sector 
C 

INCLUDE  'CONSTS. MODEL' 

C 

BO=SIGN{BO,COS(PI*TIME/22.)  ) 

Bphi=- (BO/RADIAL) * (Ws/Vsw ( POLAR) ) *SIN (POLAR) 

C 

RETURN 

END 

C 

FUNCTION  B (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

C. . .  Strength  of  the  solar  magnetic  field 
C. . .  in  a  single  sector  -  TIME- INDEPENDENT 
C 

B = SORT ( B r ( RAD I AL , POLAR , AZ I MUTH , ENERGY , T I ME ) *  *  2  + 

+  Bphi (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) **2) 

C 

C. . .  Note:  polar  component  of  B  is  identically  zero! 

C 

RETURN 

END 

C 

FUNCTION  VDSR (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

C. . .  Radial  component  of  the  drift  velocity 
C. . .  in  a  single  sector 
C 

INCLUDE  'CONSTS. MODEL' 

C 

GAMMA=1 . +ENERGY/E0 
BETA=SQRT ( 1 . -GAMMA* *  - 2 ) 

Pc =ANo*SQRT (ENERGY* (ENERGY+2 . *E0) ) 

G1=RADIAL* (Ws/Vsw (POLAR) ) *SIN (POLAR) 

G2=RADIAL* (Ws/Vsw (POLAR) )* COS (POLAR) 

C 

G=W0* SIN (POLAR) **2/ (1 . +W0*SIN(PI/2 . -POLAR) **2) 

C. . .  Note  that  G=0  if  W0=0,  i.e.,  for  a  latitude-indep .  solar  velocity 
C 

BO=SIGN(BO,COS(PI*TIME/22. )  ) 

VD0=V0* (2 . *BETA*Pc) / (3 . *ZNo*B0* (l.+Gl**2) **2) *RADIAL 
VDSR=-VD0*G2* (1 . +G* (1 . -Gl**2) ) 

C 

RETURN 

END 

C 

FUNCT I ON  VDR ( RADIAL , POLAR , AZ IMUTH , ENERGY , T I ME ) 

C. . .  Radial  component  of  the  drift  velocity 
C. . .  given  the  neutral  magnetic  sheet  SHEET 
C 

INCLUDE  ' CONSTS . MODEL ' 

C  ^  - 

GAMMA=1 . +ENERGY/E0 
BETA=SQRT (1 . -GAMMA** -2) 

Pc=ANo*SQRT (ENERGY* (ENERGY+2 . *E0) ) 

G1=RADIAL* (Ws/Vsw (POLAR) ) *SIN (POLAR) 

G2=RADIAL* (Ws/Vsw (POLAR) ) *COS (POLAR) 


+  VDSP (RADIAL, POLAR,AZIMUTH, ENERGY, TIME) **2 
+  VDSA (RADIAL , POLAR , AZIMUTH , ENERGY , TIME ) *  * 2 ) 


RETURN 

END 

FUNCTION  VD (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

Drift  speed  given  the  neutral  magnetic  sheet  SHEET 

VD=SQRT (  VDR (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) **2 
+  VDP (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) **2 
+  VDA (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) **2) 


RETURN 

END 

FUNCTION  Dt (RADIAL, POLAR , AZ IMUTH , ENERGY, TIME) 
Tangential  component  of  the  symmetric  diffusion  tensor 
kappa  (in  the  local  solar- wind  frame) 

INCLUDE  'CONSTS. MODEL' 

GAMMA=1 . +ENERGY/E0 
BETA=SQRT ( 1 . -GAMMA* *  - 2 ) 

Pc =ANo*SQRT (ENERGY* (ENERGY+2 . *E0)  ) 

RIGIDITY=Pc/ZNo 

Rel.  strength  of  local  field  to  that  at  Earth's  orbit: 
BRel=B(l. ,PI/2. , AZIMUTH, ENERGY, TIME) / 

B (RADIAL , POLAR , AZ IMUTH , ENERGY , TIME ) 

IF (RADIAL. LE. 1. )  THEN 

IF  ( RIGIDITY. LT. 1. )  THEN 
Dt=D0* (RIGIDITY) ** (2 . -EXP) *BETA 
ELSE 

Dt=0 .3*D0* (RIGIDITY) **2*BETA 
END  IF 
ELSE 

IF  (RIGIDITY. LT. 1 . )  THEN 
Dt=D0* (RIGIDITY) ** (2 . -EXP) *BETA*BRel 
ELSE 

Dt=0 .3*D0* (RIGIDITY) **2*BETA*BRel 
END  IF 
END  IF 

RETURN 

END 

FUNCTION  Dn (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

Normal  component  of  the  symmetric  diffusion  tensor 
kappa  (in  the  local  solar-wind  frame) 

INCLUDE  ' CONSTS . MODEL ' 

Dn=ETA*Dt (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

RETURN 

END 

FUNCTION  Drr (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 


noon  nonoooonn  n  o  o  onooooo  n  noonnno 


Kappa_rr  component  of  the  symmetric  diffusion  tensor  kappa 
(in  the  heliocentric  polar  spherical  coordinate  aystem) 

INCLUDE  ' CONSTS. MODEL' 

The  angle  between  the  spiral  field  line  and  the  radial  solar- wind 
direction : 

PSI=ATAN(RADIAL*Ws/Vsw( POLAR) ) 

Drr= (COS (PSD )* *2 *Dt (RADIAL, POLAR /AZIMUTH, ENERGY, TIME)  + 

+  (SIN (PSD ) * *2 *Dn (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

RETURN 

END 

FUNCTION  Drr_r (RADIAL, POLAR , AZ IMUTH , ENERGY, TIME) 

Radial  gradient  of  the  kappa_rr  component  of  the  symmetric 
diffusion  tensor  kappa 

INCLUDE  ' CONSTS . MODEL ' 

The  angle  between  the  spiral  field  line  and  the  radial  solar-wind 
direction: 

ZETA=RAD I AL* Ws / Vsw ( POLAR ) 

PSI=ATAN(ZETA) 

TERM0=2 . - (Bphi (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) / 

+  B (RADIAL, POLAR,AZIMUTH, ENERGY, TIME) ) **2 

TERM1=TERM0* (COS (PSI) **2+ETA*SIN (PSD  **2) /RADIAL 
TERM2= (Ws/ Vsw (POLAR) ) * (ETA-1. ) *SIN(2 . *PSD / (1 . +ZETA**2) 

IF (RADIAL. LE. 1. )  THEN 

Drr_r=Dt (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) * (0.  +TERM2) 

ELSE 

Drr_r=Dt (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) * (TERM1+TERM2 ) 

END  IF 

RETURN 

END 

FUNCTION  Dpp (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

Kappa_thetatheta  component  of  the  symmetric  diffusion  tensor 
kappa 

INCLUDE  ' CONSTS. MODEL' 

Note  because  kappa  is  symmetric;  Dpr=Drp=Dpa=Dap=0 . , 
these  terms  are  included  as  drift  velocity  terms! 

Dpp=Dn (RADIAL , POLAR , AZ IMUTH , ENERGY , TIME ) 

RETURN 

END 

FUNCTION  Dpp_p (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

Polar  gradient  of  the  kappa_thetethat  component  of  the  symmetric 
diffusion  tensor  kappa 


C 


INCLUDE  ' CONSTS . MODEL ' 


BO=SIGN{BO,COS(PI*TIME/22. ) ) 

TERM1=- . 5* (BO* ( Ws / Vs w (POLAR) ) /RADIAL) **2*SIN(2 . *POLAR) 

TERM2=B (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) **2 

G=W0* SIN (POLAR) **2/ ( 1 . +W0*SIN (PI/2 . -POLAR) **2) 

Note  that  G=0  if  W0=0,  i.e.,  latitude - indep .  solar-velocity ! 
TERM3=1.+2.*G 

IF (RADIAL. LE. 1. )  THEN 
Dpp_p=0 . 

ELSE 

Dpp_p=Dn (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) *TERM1*TERM3 /TERM2 
END  IF 

RETURN 

END 

FUNCTION  Daa (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 
kappa_phiphi  component  of  the  symmetric  diffusion  tensor 
Jcappa 

Note:  Daa_a  [azimuthal  gradient  of  Daa]  is  identically  zero! 
INCLUDE  ' CONSTS .MODEL' 

The  angle  between  the  spiral  field  line  and  the  radial  solar-wind 
direction : 

PSI=ATAN(RADIAL*Ws/Vsw (POLAR) ) 

Daa=COS (PSD  * *2 *Dn (RADIAL, POLAR, AZIMUTH, ENERGY, TIME)  + 

SIN (PSD  **2*Dt (RADIAL, POLAR,AZIMUTH, ENERGY, TIME) 

RETURN 

END 

FUNCTION  Dra (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 
kappa_rphi  component  of  the  symmetric  diffusion  tensor 
kappa 

INCLUDE  ' CONSTS . MODEL ' 

The  angle  between  the  spiral  field  line  and  the  radial  solar-wind 
direction: 

PSI=ATAN(RADIAL*Ws/Vsw (POLAR) ) 

Dra= (Dn (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) - 
Dt (RADIAL, POLAR,AZIMUTH, ENERGY, TIME) ) * 

COS  (PSD  *SIN(PSD 

RETURN 

END 

FUNCTION  Dra_r (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

Radial  gradient  of  the  kappa_rphi  component  of  the  symmetric 
diffusion  tensor  kappa 

INCLUDE  ' CONSTS. MODEL' 

The  angle  between  the  spiral  field  line  and  the  radial  solar-wind 
direction: 

ZETA=RADIAL*Ws/Vsw ( POLAR) 


ooono  o  o  n  o  on  o  o  non  non  n  o  non  non 


PSI=ATAN(ZETA) 

C 

TERM0=2 . - (Bphi (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) / 

+  B (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) ) **2 

TERM1=TERM0* (ETA-1. ) *SIN (2 . *PSI) / (2 . *RADIAL) 

TERM2= (ETA-1. ) * ( Ws/ Vs w (POLAR) ) *COS(2.*PSI) / (l.+ZETA**2) 

IF (RADIAL. LE. 1. )  THEN 

Dra_r=Dt (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) * (0.  +TERM2) 
ELSE 

Dra_r=Dt (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) * (TERM1+TERM2 ) 
END  IF 

RETURN 
END 

FUNCTION  Dar (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 
kappa_phir  component  of  the  symmetric  diffusion  tensor 
kappa 

INCLUDE  ' CONSTS . MODEL ' 

. .  Note  because  kappa  is  symmetric;  Dar=Dra 

Dar=Dra (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

RETURN 
END 

FUNCT I ON  Dar_r ( RADIAL , POLAR , AZ I MUTH , ENERGY , T IME ) 

Radial  gradient  of  the  kappa_phir  component  of  the  symmetric 
diffusion  tensor  kappa 

INCLUDE  ' CONSTS . MODEL ' 

Note  because  kappa  is  symmetric;  Dar_r=Dra_r 

Dar_r =Dra_r (RADIAL , POLAR , AZ IMUTH , ENERGY , TIME ) 

RETURN 
END 

FUNCTION  SHEET (RADIAL, POLAR, AZIMUTH, ENERGY, TIME) 

. .  Magnetic  current  sheet  co-rotating  with  the  Sun 

INCLUDE  ' CONSTS . MODEL ' 

POLARO=PI/2 .+TILT(TIME) *SIN (AZIMUTH+RADIAL* (Ws/Vsw ( POLAR) ) ) 

SHEET=1 . -2 . *STEP (POLAR, POLARO) 

RETURN 
END 

FUNCTION  TILT (TIME) 

Tilt  angle  of  the  Magnetic  neutral  sheet 
. .  at  the  Sun.  TILT  is  assumed  to  correlate 
with  the  11 -yr  sunspot  cycle,  i.e., 

. .  TILT  is  maximum  at  solar  sunspot  maximum  and 
TILT  is  minimum  at  solar  sunspot  minimum. 


INCLUDE  'CONSTS. MODEL' 


TILT=TILTO* (1 . + . 5*COS (2 . *PI*TIME/ll . ) ) 

Note:  With  +  in  above  instead  of  TILT  would  correspond  to 
Solar  maximum  at  t»0 ! 

RETURN 

END 

FUNCTION  STEP(X,Y) 

Approximated  Heaviside  step- function: 
for  large  INDEX 

DATA  INDEX/31/ 

STEP=.5* ( l.+TANH{ INDEX* (X-Y) ) ) • 

RETURN 

END 

FUNCTION  DIRAC (X,Y) 

Approximated  Dirac  delta- function: 
for  large  INDEX 

DATA  INDEX/31/ 

DIRAC=1. /COSH (INDEX* (X-Y) ) **2 
DIRAC=exp (- (INDEX* (X-Y) ) **2) 

RETURN 

END 


SUBROUTINE  DSS  4d (XL, XU, N1 , N2 , N3 , N4 , ND, U4D, UX4D, V) 


C.  .  .  . 

C. . .  This  is  a  straight-forward  extension  of  DSS036  to  4  independent 
C. . .  variables.  Written  by  A.  F.  Barghouty  -  July  1996. 

C .  .  . 

C. . .  SUBROUTINE  DSS_4d  COMPUTES  A  PARTIAL  DERIVATIVE  OVER  A  FOUR- 
C...  DIMENSIONAL  DOMAIN  USING  EITHER  FIVE-POINT  CENTERED  OR  FIVE- 
C...  POINT  BIASED  UPWIND  APPROXIMATIONS.  IT  IS  INTENDED  PRIMARILY 
C...  FOR  THE  NUMERICAL  METHOD  OF  LINES  (NMOL)  NUMERICAL  INTEGRATION 
C...  OF  PARTIAL  DIFFERENTIAL  EQUATIONS  (PDES)  IN  THREE  DIMENSIONS. 

C.  .  . 

C. . .  SUBROUTINE  DSS_4d  IS  CALLED  IN  ESSENTIALLY  THE  SAME  WAY  AS 
C...  SUBROUTINE  DSS036.  THE  ONLY  DIFFERENCE  IS  AN  ADDITIONAL  ARGU- 
C...  MENT,  N4,  TO  DEFINE  THE  NUMBER  OF  GRID  POINTS  IN  THE  FOURTH 
C...  DIMENSION.  THE  COMMENTS  IN  DSS036  SHOULD  THEREFORE  BE  USEFUL 
C...  IN  UNDERSTANDING  THE  OPERATION  OF  DSS_4D.  IN  PARTICULAR, 

C...  DSS036  CALLS  SUBROUTINES  DSS004  AND  DSS020  TO  IMPLEMENT  THE 
C...  FIVE-POINT  CENTERED  APPROXIMATION  AND  FIVE-POINT  BIASED  UPWIND 
C...  APPROXIMATION  OF.  THE  PARTIAL  DERIAVTIVE,  RESPECTIVELY. 

C.  .  .  ' 

C. . .  ARGUMENT  LIST 

C.  .  . 

C . . .  XL 

C.  .  . 

C.  .  . 

C . . .  XU 

C.  .  . 

C.  .  . 

C. . .  N1 

C.  .  . 

C.  .  . 

C  .  .  .  N2 

C.  .  . 

C.  .  . 

C .  .  .  N3 

C.  .  . 

C.  .  . 

C  .  .  .  N4 

C.  .  . 

C.  .  . 

C . . .  ND 

C.  .  . 

C.  .  . 

C...  U4D 

C.  .  . 

C.  .  . 

C.  .  . 

C...  UX4D 

C.  .  . 

C.  .  . 

C.  .  . 

C.  .  .  V 

C.  .  . 

C.  .  . 

C.  .  . 

C.  .  . 

C.  .  . 

C...  THE  FOLLOWING  FOUR -DIMENSIONAL  ARRAYS  CONTAIN  THE  DEPENDENT 
C...  VARIABLE  (U4D)  AND  ITS  PARTIAL  DERIVATIVE  (UX4D) 

DIMENSION  U4D (N1 , N2 , N3 , N4 ) ,  UX4D (N1 , N2 , N3 , N4 ) 


LOWER  VALUE  OF  THE  INDEPENDENT  VARIABLE  FOR  WHICH 
THE  PARTIAL  DERIVATIVE  IS  TO  BE  COMPUTED  (INPUT) 

UPPER  VALUE  OF  THE  INDEPENDENT  VARIABLE  FOR  WHICH 
THE  PARTIAL  DERIVATIVE  IS  TO  BE  COMPUTED  (INPUT) 

NUMBER  OF  GRID  POINTS  FOR  THE  FIRST  INDEPENDENT 
VARIABLE  (INPUT) 

NUMBER  OF  GRID  POINTS  FOR  THE  SECOND  INDEPENDENT 
VARIABLE  (INPUT) 

NUMBER  OF  GRID  POINTS  FOR  THE  THIRD  INDEPENDENT 
VARIABLE  (INPUT) 

NUMBER  OF  GRID  POINTS  FOR  THE  FOURTH  INDEPENDENT 
VARIABLE  (INPUT) 

NUMBER  OF  THE  INDEPENDENT  VARIABLE  FOR  WHICH  THE 
PARTIAL  DERIVATIVE  IS  TO  BE  COMPUTED  (INPUT) 

FOUR -DIMENSIONAL  ARRAY  CONTAINING  THE  DEPENDENT 
VARIABLE  WHICH  IS  TO  BE  DIFFERENTIATED  WITH  RESPECT 
TO  INDEPENDENT  VARIABLE  ND  (INPUT) 

FOUR-DIMENSIONAL  ARRAY  CONTAINING  THE  PARTIAL  DERI¬ 
VATIVE  OF  THE  DEPENDENT  VARIABLE  WITH  RESPECT  TO 
INDEPENDENT  VARIABLE  ND  (OUTPUT) 

VARIABLE  TO  SELECT  EITHER  THE  FIVE- POINT  CENTERED 
OR  FIVE-POINT  BIASED  UPWIND  APPROXIMATION  FOR  THE 
PARTIAL  DERIVATIVE.  V  EQ  0  CALLS  THE  FIVE-POINT 
CENTERED  APPROXIMATION.  V  NE  0  CALLS  THE  FIVE^POINT 
BIASED  UPWIND  APPROXIMATION  (INPUT) 
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THE  .FOLLOWING  ONE -DIMENSIONAL  ARRAYS  CONTAIN  THE  DEPENDENT 
VARIABLE  (UID)  AND  ITS  PARTIAL  DERIVATIVE  (UXID) .  IN  EACH 
CASE,  ONE  OF  THE  INDEPENDENT  VARIABLES  IS  CONSTANT  AITO  THE 
OTHER  TWO  INDEPENDENT  VARIABLES  VARY  OVER  THEIR  TOTAL  INTERVALS. 
THESE  ARRAYS  ARE  USED  FORr  TEMPORARY  STORAGE  IN  CALLING  THE 
ONE -DIMENSIONAL  ROUTINES  DSS004  AND  DSS020. 

NOTE  THAT  THE  ARRAYS  HAVE  ABSOLUTE  DIMENSIONS  AND  MAY  THERE¬ 
FORE  HAVE  TO  BE  INCREASED  IM  SIZE.  HOWEVER,  WITH  A  SIZE 
OF  51,  THE  FOUR- DIMENSIONAL  PROBLEM  COULD  HAVE  A  GRID  OF 
51  X  51  X  51  X  51  POINTS,  THEREBY  GENERATING  AN  APPROXIMATING  ODE 
SYSTEM  WITH  A  MULTIPLE  OF  51  X  51  X  51  X  51  EQUATIONS,  DEPENDING  ON 
THE  NUMBER  OF  SIMULTANEOUS  PDES .  THIS  IS  A  VERY  LARGE  ODE 
PROBLEM,  AND  THEREFORE  THE  FOLLOWING  ABSOLUTE  DIMENSIONING 
IS  CONSIDERED  ADEQUATE  FOR  MOST  PROBLEMS. 

DIMENSION  U1D(51),  UX1D(51) 

GO  TO  STATEMENT  2  IF  THE  PARTIAL  DERIVATIVE  IS  TO  BE  COMPUTED 
WITH  RESPECT  TO  THE  SECOND  INDEPENDENT  VARIABLE 
IF(ND.EQ.2)GO  TO  2 


GO  TO  STATEMENT  3  IF  THE  PARTIAL  DERIVATIVE  IS  TO  BE  COMPUTED 
WITH  RESPECT  TO  THE  THIRD  INDEPENDENT  VARIABLE 
IF(ND.EQ.3)GO  TO  3 


GO  TO  STATEMENT  4  IF  THE  PARTIAL  DERIVATIVE  IS  TO  BE  COMPUTED 
WITH  RESPECT  TO  THE  FOURTH  INDEPENDENT  VARIABLE 
IF(ND.EQ.4)GO  TO  4 


THE  PARTIAL  DERIVATIVE  IS  TO  BE  COMPUTED  WITH  RESPECT  TO  THE 
FIRST  INDEPENDENT  VARIABLE  DEFINED  OVER  AN  INTERVAL  CONSISTING 
OF  N1  GRID  POINTS.  COMPUTE  THE  PARTIAL  DERIVATIVE  AT  THE  N1  X 
N2  X  N3  X  N4  GRID  POINTS  VIA  NESTED  DO  LOOPS  09,  10,  11,  12  AND  13 
IF(Nl.EQ.l)  RETURN 
DO  10  J=1,N2 
DO  11  K=1,N3 
DO  09  L=1,N4 


TRANSFER  THE  DEPENDENT  VARIABLE  IN  THE  THREE-DIMENSIONAL  ARRAY  U3D 
TO  THE  ONE -DIMENSIONAL  ARRAY  UID  SO  THAT  SUBROUTINES  DSS004  AND 
DSS020  CAN  BE  USED  TO  CALCULATE  THE  PARTIAL  DERIVATIVE 
DO  12  1=1, N1 
UID(I)  =U4D{I,  J,K,L) 

CONTINUE 


C, 

cu, 

C, 


IF  V  EQ  0,  A  FIVE-POINT  CENTERED  APPROXIMATION  IS  USED  FOR  THE 
PARTIAL  DERIVATIVE 

IF (V.EQ. 0 . ) CALL  DSS004 (XL , XU, N1 , UID, UXID) 

IF  V  NE  0,  A  FIVE-POINT  BIASED  UPWIND  APPROXIMATION  IS  USED  FOR 
THE  PARTIAL  DERIVATIVE 

IF (V.NE. 0 . ) CALL  DSS020 (XL,XU,N1,U1D,UX1D,V) 

RETURN  THE  PARTIAL  DERIVATIVE  IN  THE  ONE -DIMENSIONAL  ARRAY  UXID 
TO  THE  FOUR -DIMENSIONAL  ARRAY  UX4D 
DO  13  1=1, N1 


UX4D(I,  J,K,L)  =UX1D(I) 

13  CONTINUE 

C.  .  . 

C. . .  THE  PARTIAL  DERIVATIVE  AT  PARTICULAR  VALUES  OF  THE  SECOND, 

C. . .  THIRD  AND  FOURTH  INDEPENDENT  VARIABLE  HAS  BEEN  CALCULATED. 

C. . .  REPEAT  THE  CALCULATION  FOR  THE  OTHER  VALUES  OF  THE  SECOND, 

C...  THIRD,  AND  FOURTH  INDEPENDENT  VARIABLES 
09  CONTINUE 

11  CONTINUE 

10  CONTINUE 

C.  .  . 

C . . .  THE  PARTIAL  DERIVATIVE  HAS  BEEN  CALCULATED  OVER  THE  ENTIRE  N1  X 
C. . ,  N2  X  N3  X  N4  GRID.  THEREFORE  RETURN  TO  THE  CALLING  PROGRAM  WITH 
C...  THE  PARTIAL  DERIVATIVE  IN  THE  FOUR -DIMENSIONAL  ARRAY  UX4D 
RETURN 

C.  .  . 

Q  ★tIt***********************************************************'*'*'*** 

c. . . 

C...  THE  PARTIAL  DERIVATIVE  IS  TO  BE  COMPUTED  WITH  RESPECT  TO  THE 
C.  ’.  .  SECOND  INDEPENDENT  VARIABLE  DEFINED  OVER  AN  INTERVAL  CONSISTING 
C . . .  OF  N2  GRID  POINTS .  COMPUTE  THE  PARTIAL  DERIVATIVE  AT  THE  N1  X 
C. . .  N2  X  N3  X  N4  GRID  POINTS  VIA  NESTED  DO  LOOPS  19,  20,  21,  22  AND  23 
2  IF{N2.EQ.l)  RETURN 

DO  20  1=1, N1 
DO  21  K=1,N3 
DO  19  L=1,N4 

C.  .  . 

C...  TRANSFER  THE  DEPENDENT  VARIABLE  IN  THE  FOUR -DIMENSIONAL  ARRAY  U4D 
C. . .  TO  THE  ONE-DIMENSIONAL  ARRAY  UID  SO  THAT  SUBROUTINES  DSS004  AND 
C. . .  DSS020  CAN  BE  USED  TO  CALCULATE  THE  PARTIAL  DERIVATIVE 
DO  22  J=1,N2 
UID(J) =U4D(I, J,K,L) 

22  CONTINUE 
C.  .  . 

C...  IF  V  EQ  0,  A  FIVE-POINT  CENTERED  APPROXIMATION  IS  USED  FOR  THE 
C. . .  PARTIAL  DERIVATIVE 

IF (V.EQ. 0 . ) CALL  DSS004 (XL, XU, N2 , UID, UXID) 

C.  .  . 

C...  IF  V  NE  0,  A  FIVE-POINT  BIASED  UPWIND  APPROXIMATION  IS  USED  FOR 
C...  THE  PARTIAL  DERIVATIVE 

IF (V.NE. 0 . ) CALL  DSS020 (XL, XU, N2 , UID, UXID, V) 

C.  .  . 

C...  RETURN  THE  PARTIAL  DERIVATIVE  IN  THE  ONE -DIMENSIONAL  ARRAY  UXID 
C...  TO  THE  FOUR -DIMENSIONAL  ARRAY  UX4D 
DO  23  J=1,N2 
UX4D(I,  J,K,L)  =UX1D(J) 

23  CONTINUE 
C.  .  . 

C. . .  THE  PARTIAL  DERIVATIVE  AT  PARTICULAR  VALUES  OF  THE  FIRST, 

C...  THIRD,  AND  FOURTH  INDEPENDENT  VARIABLE  HAS  BEEN  CALCULATED. 

C. . .  REPEAT  THE  CALCULATION  FOR  THE  OTHER  VALUES  OF  THE  FIRST, 

C...  THIRD,  AND  FOURTH  INDEPENDENT  VARIABLES 

19  CONTINUE 

21  CONTINUE 

20  CONTINUE 

C.  .  . 

C...  THE  PARTIAL  DERIVATIVE  HAS  BEEN  CALCULATED  OVER  THE  ENTIRE  N1  X 
C. . .  N2  X  N3  X  N4  GRID.  THEREFORE  RETURN  TO  THE  CALLING  PROGRAM  WITH 
C...  THE  PARTIAL  DERIVATIVE  IN  THE  FOUR -DIMENSIONAL  ARRAY  UX4D 
RETURN 


c 


THE  PARTIAL  DERIVATIVE  IS  TO  BE  COMPUTED  WITH  RESPECT  TO  THE 
THIRD  INDEPENDENT  VARIABLE  DEFINED  OVER  AN  INTERVAL  CONSISTING 
OF  N3  GRID  POINTS.  COMPUTE  THE  PARTIAL  DERIVATIVE  AT  THE  N1  X 
N2  X  N3  X  N4  GRID  POINTS  VIA  NESTED  DO  LOOPS  29,  30,  31,  32  AND  33 
IF(N3.EQ.l)  RETURN 
DO  30  1=1, N1 
DO  31  J=1,N2 
DO  29  L=1,N4 


TRANSFER  THE  DEPENDENT  VARIABLE  IN  THE  FOUR-DIMENSIONAL  ARRAY  U4D 
TO  THE  ONE-DIMENSIONAL  ARRAY  UID  SO  THAT  SUBROUTINES  DSS004  AND 
DSS020  CAN  BE  USED  TO  CALCULATE  THE  PARTIAL  DERIVATIVE 
DO  32  K=1,N3 
UID(K)  =U4D(I,  J,K,L) 

CONTINUE 

IF  V  EQ  0,  A  FIVE-POINT  CENTERED  APPROXIMATION  IS  USED  FOR  THE 
PARTIAL  DERIVATIVE 

IF (V.EQ. 0 . ) CALL  DSS004 (XL, XU, N3 , UID, UXID) 

IF  V  NE  0,  A  FIVE-POINT  BIASED  UPWIND  APPROXIMATION  IS  USED  FOR 
THE  PARTIAL  DERIVATIVE 

IF (V.NE. 0 . ) CALL  DSS020 (XL, XU, N3 , UID, UXID, V) 

RETURN  THE  PARTIAL  DERIVATIVE  IN  THE  ONE -DIMENSIONAL  ARRAY  UXID 
TO  THE  FOUR -DIMENSIONAL  ARRAY  UX4D 
DO  33  K=1,N3 
UX4D(I, J,K,L) =UX1D(K) 

CONTINUE 

THE  PARTIAL  DERIVATIVE  AT  PARTICULAR  VALUES  OF  THE  FIRST, 

SECOND,  AND  FOURTH  INDEPENDENT  VARIABLE  HAS  BEEN  CALCULATED. 

REPEAT  THE  CALCULATION  FOR  THE  OTHER  VALUES  OF  THE  FIRST, 

SECOND  AND  FOURTH  INDEPENDENT  VARIABLES 

CONTINUE 

CONTINUE 

CONTINUE 

THE  PARTIAL  DERIVATIVE  HAS  BEEN  CALCULATED  OVER  THE  ENTIRE  N1  X 
N2  X  N3  X  N4  GRID.  THEREFORE  RETURN  TO  THE  CALLING  PROGRAM  WITH  THE 
PARTIAL  DERIVATIVE  IN  THE  FOUR -DIMENSIONAL  ARRAY  UX4D 
RETURN 

***★**★★★**★★★★★★★★★★★★***★*★***★*★**★***★★★★**★★*★****'★*★★■*’■*■★★*★★ 

THE  PARTIAL  DERIVATIVE  IS  TO  BE  COMPUTED  WITH  RESPECT  TO  THE 
FOURTH  INDEPENDENT  VARIABLE  DEFINED  OVER  AN  INTERVAL  CONSISTING 
OF  N4  GRID  POINTS.  COMPUTE  THE  PARTIAL  DERIVATIVE  AT  THE  N1  X 
N2  X  N3  X  N4  GRID  POINTS  VIA  NESTED  DO  LOOPS  40,  41,  42,  43  AND  44 
IF(N4.EQ.l)  RETURN 
DO  40  1=1, N1 
DO  41  J=1,N2 
DO  42  K=1,N3 


TRANSFER  THE  DEPENDENT  VARIABLE  IN  THE  FOUR -DIMENSIONAL  ARRAY  U4D 
TO  THE  ONE -DIMENSIONAL  ARRAY  UID  SO  THAT  SUBROUTINES  DSS004  AND 
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DSS020  CAN  BE  USED  TO  CALCULATE  THE  PARTIAL  DERIVATIVE 
DO  43  L=1,N4 
UID(L) =U4D(I, J,K,L) 

CONTINUE 

IF  V  EQ  0,  A  FIVE-POINT  CENTERED  APPROXIMATION  IS  USED  FOR  THE 
PARTIAL  DERIVATIVE 

IF (V.EQ. 0 . ) CALL  DSS004 (XL, XU, N4 , UID, UXID) 

IF  V  NE  0,  A  FIVE-POINT  BIASED  UPWIND  APPROXIMATION  IS  USED  FOR 
THE  PARTIAL  DERIVATIVE 

IF(V.NE. 0 . ) CALL  DSS020 (XL, XU, N4 , UID, UXID, V) 

RETURN  THE  PARTIAL  DERIVATIVE  IN  THE  ONE -DIMENSIONAL  ARRAY  UXID 
TO  THE  FOUR -DIMENSIONAL  ARRAY  UX4D 
DO  44  L=1,N4 
UX4D(I, J,K,L) =UX1D(L) 

CONTINUE 

THE  PARTIAL  DERIVATIVE  AT  PARTICULAR  VALUES  OF  THE  FIRST, 

SECOND,  AND  THIRD  INDEPENDENT  VARIABLE  HAS  BEEN  CALCULATED. 

REPEAT  THE  CALCULATION  FOR  THE  OTHER  VALUES  OF  THE  FIRST, 

SECOND  AND  THIRD  INDEPENDENT  VARIABLES 

CONTINUE 

CONTINUE 

CONTINUE 

THE  PARTIAL  DERIVATIVE  HAS  BEEN  CALCULATED  OVER  THE  ENTIRE  N1  X 

N2  X  N3  X  N4  GRID.  THEREFORE  RETURN  TO  THE  CALLING  PROGRAM  WITH  THE 

PARTIAL  DERIVATIVE  IN  THE  FOUR -DIMENSIONAL  ARRAY  UX4D 

RETURN 

END 

SUBROUTINE  DSS004 (XL, XU, N, U, UX) 

SUBROUTINE  DSS004  COMPUTES  THE  FIRST  DERIVATIVE,  U  ,  OF  A 

X 

VARIABLE  U  OVER  THE  SPATIAL  DOMAIN  XL  LE  X  LE  XU  FROM  CLASSICAL 
FIVE-POINT,  FOURTH-ORDER  FINITE  DIFFERENCE  APPROXIMATIONS 

ARGUMENT  LIST 

XL  LOWER  BOUNDARY  VALUE  OF  X  (INPUT) 

XU  UPPER  BOUNDARY  VALUE  OF  X  (INPUT) 

N  NUMBER  OF  GRID  POINTS  IN  THE  X  DOMAIN  INCLUDING  THE 

BOUNDARY  POINTS  (INPUT) 

U  ONE-DIMENSIONAL  ARRAY  CONTAINING  THE  VALUES  OF  U  AT 

THE  N  GRID  POINT  POINTS  FOR  WHICH  THE  DERIVATIVE  IS 
TO  BE  COMPUTED  (INPUT) 

UX  ONE-DIMENSIONAL  ARRAY  CONTAINING  THE  NUMERICAL  ' 

VALUES  OF  THE  DERIVATIVES  OF  U  AT  THE  N  GRID  POINTS  - 
(OUTPUT) 

THE  MATHEMATICAL  DETAILS  OF  THE  FOLLOWING  TAYLOR  SERIES  (OR 
POLYNOMIALS)  ARE  GIVEN  IN  SUBROUTINE  DSS002 . 


|n(  pn  _/Jio  Jo  r  fio  Jot  ho  Jo  jor.J^o  Joc  J~)  o .  _  J  o  Jo  r. 


FIVE- POINT  FORMULAS 


I 

I 


C 

c., 

C 


(1)  LEFT  END,  POINT 

1  =  1 

2 

3 

4 

A(U2  =  U1  +  U1  (  DX) 

+ 

U1  (  DX) 

+ 

U1  (  DX) 

+ 

U1 

(  DX) 

X  IF 

2X  2F 

3X  3F 

4X 

4F 

5 

6 

7 

+  U1  (  DX) 

+ 

U1  J  DX) 

+ 

U1  (  DX) 

+ 

.  .  .  ) 

5X  5F 

6X  6F 

7X  7F 

2 

3 

4 

B(U3  =  U1  +  U1  (2DX) 

+ 

U1  {2DX) 

+ 

U1  (2DX) 

+ 

U1 

{2DX) 

X  IF 

2X  2F 

3X  3F 

4X 

4F 

5 

6 

7 

+  U1  (2DX) 

+ 

U1  (2DX) 

+ 

U1  (2DX) 

+ 

.  .  .  ) 

5X  5F 

6X  6F 

7X  7F 

2 

3 

4 

C(U4  =  U1  +  U1  (3DX) 

+ 

U1  (3DX) 

+ 

U1  (3DX) 

+ 

U1 

(3DX) 

X  IF 

2X  2F 

3X  3F 

4X 

4F 

5 

6 

7 

+  U1  (3DX) 

+ 

U1  (3DX) 

+ 

U1  {3DX) 

+ 

.  .  .  ) 

5X  5F 

6X  6F 

7X  7F 

2 

3 

4 

D(U5  =  U1  +  U1  (4DX) 

+ 

U1  (4DX) 

+ 

U1  (4DX) 

+ 

U1 

(4DX) 

X  IF 

2X  2F 

3X  3F 

4X 

4F 

5 

6 

7 

+  U1  (4DX) 

+ 

U1  (4DX) 

+ 

U1  (4DX) 

+ 

.  .  .) 

5X  5F 

6X  6F 

7X  7F 

CONSTANTS  A,  B,  C  AND  D  ARE  SELECTED  SO  THAT  THE  COEFFICIENTS 
OF  THE  U1  TERMS  SUM  TO  ONE  AND  THE  COEFFICIENTS  OF  THE  U1  , 


X 


U1 

3X 

AND 

U1 

TERMS  SUM 
4X 

TO 

A 

+ 

2B 

+ 

3C 

+ 

4D 

=  1 

A 

+ 

4B 

+ 

9C 

+ 

16D 

=  0 

A 

+ 

8B 

+ 

27C 

+ 

64D 

=  0 

A 

+ 

16B 

+ 

81C 

+ 

256D 

=  0 

2X 


SIMULTANEOUS  SOLUTION  FOR  A,  B,  C  AND  D  FOLLOWED  BY  THE  SOLU¬ 
TION  OF  THE  PRECEDING  TAYLOR  SERIES,  TRUNCATED  AFTER  THE  U 

4X 

TERMS,  FOR  U1  GIVES  THE  FOLLOWING  FIVE-POINT  APPROXIMATION 
X 


4 

U1  =  (1/12DX) (-25U1  +  48U2  -  36U3  +  16U4  -  3U5)  +  0(DX  )  (1) 

X 


C 


(2)  INTERIOR  POINT,  1=2 


. 

2 

3 

4 

A(U1  =  U2  +  U2  (-DX) 

+ 

U2  (-DX) 

+ 

U2  (-DX) 

+ 

U2  (-DX) 

X  IF 

2X  2F 

3X  3F 

4X  4F 

5 

6 

7 

+  U2  (-DX) 

+ 

U2  ( -DX) 

+ 

U2  (-DX) 

+ 

.  .  .) 

5X  5F 

6X  6F 

7X  7F 

2 

3 

4 

B(U3  =  U2  +  U2  (  DX) 

+ 

U2  (  DX) 

+ 

U2  (  DX) 

+ 

U2  (  DX) 

X  IF 

2X  2F  ■ 

3X  3F 

4X  4F 

5 

6 

7 

+  U2  (  DX) 

+ 

U2  (  DX) 

+ 

U2  (  DX) 

+ 

.  .  .) 

5X  5F 

6X  6F 

7X  7F 

2 

3 

4 

C(U4  =  U2  +  U2  (2DX) 

+ 

U2  (2DX) 

+ 

U2  (2DX) 

+ 

U2  (2DX) 

X  IF 

2X  2F 

3X  3F 

4X  4F 

5 

6 

7 

+  U2  (2DX) 

+ 

U2  (2DX) 

+ 

U2  (2DX) 

+ 

.  .  .  ) 

5X  5F 

6X  6F 

7X  7F 

2 

3 

4 

D(U5  =  U2  +  U2  (3DX) 

+ 

U2  (3DX) 

+ 

U2  (3DX) 

+ 

U2  (3DX) 

X  IF 

2X  2F 

3X  3F 

4X  4F 

5 

6 

7 

+  U2  (3DX) 

+ 

U2  (3DX) 

+ 

U2  (3DX) 

+ 

.  .  .) 

5X  5F 

6X  6F 

7X  7F 

-A  +  B  +  2C  +  3D  =  1 


A+  B  +  4C+  9D  =  0 
-A  +  B  +  8C  +  27D  =  0 
A  +  B  +  16C  +  81D  =  0 


SIMULTANEOUS  SOLUTION  FOR  A,  B,  C  AND  D  FOLLOWED  BY  THE  SOLU¬ 
TION  OF  THE  PRECEDING  TAYLOR  SERIES,  TRUNCATED  AFTER  THE  U 

4X 

TERMS,  FOR  U1  GIVES  THE  FOLLOWING  FIVE -POINT  APPROXIMATION 
X 

4 

U2  =  (1/12DX) (-3U1  -  10U2  +  18U3  -  6U4  +  U5)  +  0(DX  )  (2) 

X 

(3)  INTERIOR  POINT  I,  I  NE  2,  N-1 

2  3 

A(UI-2  =  UI  +  UI  (-2DX)  +  UI  (-2DX)  +  UI  (-2DX) 

X  IF  2X  2F  3X  3F 


4  5  6 

+  UI  (-2DX)  +  UI  (-2DX)  +  UI  (-2DX)  + 

4X  4F  5X  5F  6X  6F 


1 

2 

3 

c. .  . 

B  (UI-1 

=  UI  +  UI  (  -DX)  + 

UI 

(  -DX) 

+  UI  ( 

-DX) 

-«=) 

X  IF 

2X 

2F 

3X 

3F 

c. .  . 

4 

5 

6 

■  ■ 

+  UI  (  -DX)  + 

UI 

(  -DX) 

+  UI  ( 

-DX)  + 

.  .  .  ) 

' 

4X  4F 

5X 

5F 

6X 

6F 

c. .  . 

2 

3 

.  . 

C(UI  +  1 

=  UI  +  UI  (  DX)  + 

TJI 

(  DX) 

+  UI  ( 

DX) 

' .  .  . 

X  IF 

2X 

2F 

3X 

3F 

C.  .  . 

.  . 

4 

5 

• 

6 

+  UI  {  DX)  + 

UI 

(  DX) 

+  UI  ( 

DX)  + 

.  .  .  ) 

^  . 

4X  4F 

5X 

5F 

6X 

6F 

C.  .  . 

2 

3 

•  ^  • 

D(UI+2 

=  UI  +  UI  (  2DX)  + 

UI 

(  2DX) 

+  UI  ( 

2DX) 

c. .  . 

X  IF 

2X 

2F 

3X 

3F 

■C.  .  . 

4 

5 

6 

.  .  . 

+  UI  (  2DX)  + 

UI 

{  2DX) 

+  UI  ( 

2DX)  + 

.  .  .  ) 

£■  ■  • 

4X  4F 

5X 

5F 

6X 

6F 

-2A  - 

B  +  C  +  2D  =  1 

c. . . 

4A  + 

B  +  C  +  4D  =  0 

c. .  . 

-8A  - 

B  +  C  +  8D  =  0 

.  . 

16A  + 

B  +  C  +  16D  =  0 

C.  .  . 

SIMULTANEOUS  SOLUTION  FOR  A, 

B, 

C  AND  D 

FOLLOWED  BY  THE 

SOLU 

TION  OF  THE  PRECEDING  TAYLOR  SERIES,  TRUNCATED  AFTER  THE 

U 

4X 

c. .  . 

TERMS , 

FOR  UI  GIVES  THE  FOLLOWING  FIVE 

-POINT  APPROXIMATION 

•C'.  .  . 

X 

4 

UI  = 

(1/12DX) (UI-2  -  8UI-1 

+  GUI  +  8UI+1  -  UI+2)  +  0{DX 

) 

c. .  . 

X 

(4) 

INTERIOR  POINT,  I  = 

N-1 

c. .  . 

A 

2 

.  . 

A(UN-4 

=  UN-1  +  UN-1  (-3DX) 

+ 

UN-1  ( 

-3DX)  + 

UN-1  (- 

3DX) 

c .  .  . 

X  IF 

2X 

2F 

3X 

3F 

.  . 

4 

5 

6 

w .  .  . 

+  UN-1  (-3DX)  +  UN- 

1  ( 

-3DX) 

+  UN-1 

(-3DX)  + 

•  .  • 

c. .  . 

4X  4F 

5X 

5F 

6X 

6F 

2 

c. .  . 

B  (UN-3 

=  UN-1  +  UN-1  (-2DX) 

+ 

UN-1  ( 

-2DX)  + 

UN-1  (- 

2DX) 

A. 

X  IF 

2X 

2F 

3X 

3F 

^ .  .  . 

4 

5 

6 

£•  •  • 

+  UN-1  (-2DX)  +  UN- 

1  ( 

-2DX) 

+  UN-1 

(-2DX)  + 

•  •  • 

4X  4F 

5X 

5F 

6X 

6F 

c. . . 

2 

c. .  . 

C(UN- 

2  =  UN- 

1  + 

UN-1  (  -DX)  + 

UN-1  (- 

-X)  + 

UN-1  ( 

c. .  . 

. 

X  IF 

2X 

2F 

3X 

c. .  . 
c. .  . 

4 

5 

6 

c. .  . 

+  UN- 

1 

( 

-DX)  +  UN-1 

(  -DX)  + 

UN-1 

(  -DX)  + 

c. .  . 
c. .  . 

4X 

4F  5X 

5F 

6X 

6F 

c. .  . 

2 

c. .  . 

D  (UN 

=  UN- 

1  + 

UN-1  (  DX)  + 

UN-1  ( 

DX)  + 

UN-1  ( 

c. .  . 

X  IF  - 

2X 

2F 

3X 

c. .  . 
c. .  . 

4 

5 

6 

c. .  . 

+  UN- 

1 

( 

DX)  +  UN-1 

(  DX)  + 

UN-1 

(  DX)  + 

c. .  . 
c. .  . 

4X 

4F  5X 

5F 

6X 

6F 

c. .  . 

-3A  - 

2B  - 

C 

+ 

D  =  1 

c. .  . 

c. .  . 

9A  + 

4B  + 

c 

+ 

o 

II 

Q 

c. .  . 

c. .  . 

-27A  - 

8B  - 

c 

+ 

o 

II 

Q 

c. .  . 
c. .  . 

81A  + 

16B  + 

c 

+ 

o 

II 

Q 

-DX) 


C.  .  . 
C.  .  . 
C.  .  . 
C.  .  . 
C.  .  . 
C.  .  . 
C.  .  . 

c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  , 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 


SIMULTANEOUS  SOLUTION  FOR  A,  B,  C  AND  D  FOLLOWED  BY  THE  SOLU¬ 
TION  OF  THE  PRECEDING  TAYLOR  SERIES,  TRUNCATED  AFTER  THE  U 

4X 

TERMS,  FOR  U1  GIVES  THE  FOLLOWING  FIVE -POINT  APPROXIMATION 
X 

4 

UN-1  =  (1/12DX) (-UN-4  +  6UN-3  -  18UN-2  +  lOUN-1  +  3UN)  +  0(DX  ) 
X 

(4) 

(5)  RIGHT  END,  POINT  I  =  N 


A(UN-4  =  UN  +  UN  (-4DX)  +  UN  (-4DX)  +  UN  (-4DX) 

X  IF  2X  2F  3X  3F 


+  UN  (-4DX)  +  UN  (-4DX)  +  UN  (-4DX)  + 

4X  4F  5X  5F  6X  6F 


B(UN-3  =  UN  +  UN  {-3DX)  +  UN  (-3DX)  +  UN  (-3DX) 

X  IF  2X  2F  3X  3F 


UN  (-3DX)  +  UN  (-3DX)  +  UN  (-3DX)  + 

4X  4F  5X  5F  6X  6F 


C(UN-2  =  UN  +  UN  (-2DX)  +  UN  (-2DX)  +  UN  (-2DX) 

X  IF  2X  2F  3X  3F 


+  UN  (-2DX)  +  UN  (-2DX)  +  UN  (-2DX)  + 

4X  4F  5X  5F  6X  6F 


.  ) 


2 


3 


po  |o  r  pn  jor  ho  in  {o  r  p  o  In  f  h  n  loi  ho 


D(UN-1 

=  UN 

+  UN  ( 

-DX) 

- 

X 

IF 

4 

+ 

UN  ( 

-DX) 

4X 

4F 

-4A  - 

3B  - 

2C  - 

D  = 

16A  + 

9B  + 

4C  + 

D  = 

-64A  - 

27B  - 

8C  - 

D  = 

256A  + 

81B  + 

16C  + 

D  = 

+  UN  ( 

-DX) 

+  UN  ( 

-DX) 

2X 

2F 

3X 

3F 

5 

+  UN  ( 

-DX) 

+  UN  ( 

-DX) 

5X 

5F 

6X 

6F 

1 

0 

0 

0 


SIMULTANEOUS  SOLUTION  FOR  A,  B,  C  AND  D  FOLLOWED  BY  THE  SOLU¬ 
TION  OF  THE  PRECEDING  TAYLOR  SERIES,  TRUNCATED  AFTER  THE  U 

4X 

TERMS,  FOR  U1  GIVES  THE  FOLLOWING  FIVE-POINT  APPROXIMATION 
X 


4 

UN  =  (1/12DX) (3UN-4  -  16UN-3  +  36UN-2  -  48UN-1  +  25UN)  +  0(DX  ) 


X 

(5) 


THE  WEIGHTING  COEFFICIENTS  FOR  EQUATIONS  (1)  TO  (5)  CAN  BE 
SUMMARIZED  AS 


-25 

48 

-36 

16 

-3 

-3 

-10 

18 

-6 

1 

1 

-8 

0 

8 

-1 

-1 

6 

-18 

10 

3 

3 

-16 

36 

-48 

25 

WHICH  ARE  THE  COEFFICIENTS  REPORTED  BY  BICKLEY  FOR  N  =  4,  M  = 

1,  P=0,  1,  2,  3,  4  (BICKLEY,  W.  G. ,  FORMULAE  FOR  NUMERICAL 
DIFFERENTIATION,  MATH.  GAZ . ,  VOL.  25,  1941.  NOTE  -  THE  BICKLEY 
COEFFICIENTS  HAVE  BEEN  DIVIDED  BY  A  COMMON  FACTOR  OF  TWO) . 

EQUATIONS  (1)  TO  (5)  CAN  NOW  BE  PROGRAMMED  TO  GENERATE  THE 
DERIVATIVE  U  (X)  OF  FUNCTION  U(X)  (ARGUMENTS  U  AND  UX  OF  SUB- 
X 

ROUTINE  DSS004  RESPECTIVELY) . 


C 


c , 


DIMENSION  U(N),UX{N) 

COMPUTE  THE  SPATIAL  INCREMENT 
DX= (XU-XL) /FLOAT (N-1) 

R4FDX=1./ (12.*DX) 

NM2=N-2 

EQUATION  (1)  (NOTE  -  THE  RHS  OF  EQUATIONS  (1)',  (2),  (3),  (4) 

AND  (5)  HAVE  BEEN  FORMATTED  SO  THAT  THE  NUMERICAL  WEIGHTING 
COEFFICIENTS  CAN  BE  MORE  EASILY  ASSOCIATED  WITH  THE  BICKLEY 
MATRIX  ABOVE) 

UX(  1)=R4FDX* 


on  nooo  ooononooonoono  n  oo  non  no  on 


1(  -25.*U(  1)  +48.*U(  2)  -36.*U(  3)  +16.*U(  4)  -3.*U{  5)) 

EQUATION  (2) 

UX(  2)=R4FDX* 

1(  -3.*U{  1)  -10.*U(  2)  +18.*U(  3)  -6.*U(  4)  +1.*U(  5)) 

EQUATION  (3) 

DO  1  1=3, NM2 
UX{  I)=R4FDX* 

1(  +l.*U(I-2)  -8.*U(I-1)  +D.*U(  I)  +8.*U{I+1)  -l.*U(I+2)) 

CONTINUE 

EQUATION  (4) 

UX (N-1) =R4FDX* 

1(  -l.*U(N-4)  +6.*U(N-3)  -18.*U(N-2)  +10.*U(N-1)  +3.*U(  N) ) 

EQUATION  (5) 

UX (  N) =R4FDX* 

1(  3.*U(N-4)  -16.*U(N-3)  +36.*U(N-2)  -48.*U{N-1)  +25.*U(  N) ) 

RETURN 
END 

SUBROUTINE  DSS020 (XL, XU, N, U, UX, V) 

SUBROUTINE  DSS020  IS  AN  APPLICATION  OF  FOURTH-ORDER  DIRECTIONAL 
DIFFERENCING  IN  THE  NUMERICAL  METHOD  OF  LINES.  IT  IS  INTENDED 
SPECIFICALLY  FOR  THE  ANALYSIS  OF  CONVECTIVE  SYSTEMS  MODELLED  BY 
FIRST-ORDER  HYPERBOLIC  PARTIAL  DIFFERENTIAL  EQUATIONS  AS  DIS¬ 
CUSSED  IN  SUBROUTINE  DSS012.  THE  COEFFICIENTS  OF  THE  FINITE 
DIFFERENCE  APPROXIMATIONS  USED  HEREIN  ARE  TAKEN  FROM  BICKLEY,  W. 
G.,  FORMULAE  FOR  NUMERICAL  DIFFERENTIATION,  THE  MATHEMATICAL 
GAZETTE,  PP.  19-27,  1941,  N  =  4,  M  =  1,  P  =  0,  1,  2,  3,  4.  THE 
IMPLEMENTATION  IS  THE  **FIVE- POINT  BIASED  UPWIND  FORMULA**  OF 
M.  B.  CARVER  AND  H.  W.  HINDS,  THE  METHOD  OF  LINES  AND  THE 
ADVECTION  EQUATION,  SIMULATION,  VOL.  31,  NO.  2,  PP.  59-69, 

AUGUST,  1978 

DIMENSION  U(N),UX(N) 

COMPUTE  THE  COMMON  FACTOR  FOR  EACH  FINITE  DIFFERENCE  APPROXIMATION 
CONTAINING  THE  SPATIAL  INCREMENT,  THEN  SELECT  THE  FINITE  DIFFER¬ 
ENCE  APPROXIMATION  DEPENDING  ON  THE  SIGN  OF  V  (SIXTH  ARGUMENT) . 

DX= (XU-XL) /FLOAT (N-1) 

R4FDX=1./ (12.*DX) 

IF (V.LT. 0 . ) GO  TO  10 

(1)  FINITE  DIFFERENCE  APPROXIMATION  FOR  POSITIVE  V 
UX(  1)=R4FDX* 


1(  -25.*U(  1) 

UX(  2)=R4FDX* 

+48.*U(  2) 

-36.*U(  3) 

+16 . *U( 

4) 

-3 . *U( 

5)  ) 

1(  -3.*U(  1) 

UX(  3)=R4FDX* 

-10.*U(  2) 

+18.*U(  3) 

-6 . *U( 

4) 

+1 . *U( 

5)  ) 

1(  +1.*U(  1) 

NM1=N-1 

DO  1  1=4, NMl 
UX(  I)=R4FDX* 

-8.*U(  2) 

+0.*U(  3) 

+  8  . *U( 

4) 

-l.*U( 

5)  ) 

1(  -l.*U(I-3) 

CONTINUE 

UX(  N)=R4FDX* 

+6 . *U(I-2) 

-18 . *U(I-1) 

+10 . *U( 

I) 

+3 . *U(I+1) ) 

1(  3.*U(N-4) 

-16 . *U(N-3) 

+36 . *U(N-2) 

-48 . *U(N- 

-1) 

+25 . *U( 

N)  ) 

RETURN 


(2)  FINITE  DIFFERENCE 
UX(  1)=R4FDX* 

1(  -25.*U(  1)  +48.*U(  2) 

NM3=N-3 
DO  2  1=2, NM3 
UX(  I)=R4FDX* 

1(  -3.*U(I-1)  -10.*U(  I) 

CONTINUE 
UX(N-2) =R4FDX* 

1(  +l.*U(N-4)  -8.*U(N-3) 

UX(N-l) =R4FDX* 

1(  -l.*U(N-4)  +6.*U(N-3) 

UX(  N) =R4FDX* 

1(  3.*U(N-4)  -16.*U(N-3) 

RETURN 
END 


APPROXIMATION  FOR  NEGATIVE  V 


-36.*U(  3) 

+18.*U(I+1) 

+0.*U(N-2) 
-18 . *U(N-2) 
+36 . *U{N-2) 


+16.*U{  4) 

-6. *U(I+2) 

+8 . *U(N-1) 
+10 . *U(N-1) 
-48, *U(N-1) 


-3.*U(  5)) 

+l.*U(I+3) ) 

-l.*U(  N)  ) 
+3.*U(  N) ) 
+25.*U(  N)  ) 


SUBROUTINE  LSODES  (F,  NEQ,  Y,  T,  TOUT,  ITOL,  RTOL,  ATOL,  ITASK, 

1  ISTATE,  lOPT,  RWORK,  LRW,  IWORK,  LIW,  JAC,  MF) 

EXTERNAL  F,  JAC 

INTEGER  NEQ,  ITOL,  ITASK,  ISTATE,  lOPT,  LRW,  IWORK,  LIW,  MF 
REAL  Y,  T,  TOUT,  RTOL,  ATOL,  RWORK 

DIMENSION  NEQ (1) ,  Y{1),  RTOL(l),  ATOL(l),  RWORK(LRW),  IWORK(LIW) 

C - 

C  THIS  IS  THE  MAY  2,  1983  VERSION  OF 

C  LSODES . .  LIVERMORE  SOLVER  FOR  ORDINARY  DIFFERENTIAL  EQUATIONS 
C  WITH  GENERAL  SPARSE  JACOBIAN  MATRICES. 

C  THIS  VERSION  IS  IN  SINGLE  PRECISION. 

C 

C  LSODES  SOLVES  THE  INITIAL  VALUE  PROBLEM  FOR  STIFF  OR  NONSTIFF 
C  SYSTEMS  OF  FIRST  ORDER  ODE-S, 

C  DY/DT  =  F(T,Y)  ,  OR,  IN  COMPONENT  FORM, 

C  DY(I)/DT  =  F{I)  =  F(I,T,Y(1)  ,Y(2)  ,  .  .  .  ,Y{NEQ)  )  (I  =  1 _ ,NEQ). 

C  LSODES  IS  A  VARIANT  OF  THE  LSODE  PACKAGE,  AND  IS  INTENDED  FOR 
C  PROBLEMS  IN  WHICH  THE  JACOBIAN  MATRIX  DF/DY  HAS  AN  ARBITRARY 
C  SPARSE  STRUCTURE  (WHEN  THE  PROBLEM  IS  STIFF) . 

C 

C  AUTHORS..  ALAN  C.  HINDMARSH, 

C  MATHEMATICS  AND  STATISTICS  DIVISION,  L-316 

C  LAWRENCE  LIVERMORE  NATIONAL  LABORATORY 

C  LIVERMORE,  CA  94550. 

C 

C  AND  ANDREW  H.  SHERMAN 

C  EXXON  PRODUCTION  RESEARCH  CO. 

C  P.  0.  BOX  2189 

C  HOUSTON,  TX  77001 

C - - - 

C  REFERENCES . . 

C  1.  ALAN  C.  HINDMARSH,  LSODE  AND  LSODI,  TWO  NEW  INITIAL  VALUE 
C  ORDINARY  DIFFERENTIAL  EQUATION  SOLVERS, 

C  ACM- SIGNUM  NEWSLETTER,  VOL.  15,  NO.  4  (1980),  PP .  10-11. 

C 

C  2.  S.  C.  EISENSTAT,  M.  C.  GURSKY,  M.  H.  SCHULTZ,  AND  A.  H.  SHERMAN, 

C  YALE  SPARSE  MATRIX  PACKAGE..  I.  THE  SYMMETRIC  CODES, 

C  RESEARCH  REPORT  NO.  112,  DEPT.  OF  COMPUTER  SCIENCES,  YALE 

C  UNIVERSITY,  1977. 

C 

C  3.  S.  C.  EISENSTAT,  M.  C.  GURSKY,  M.  H.  SCHULTZ,  AND  A.  H.  SHERMAN, 

C  YALE  SPARSE  MATRIX  PACKAGE..  II.  THE  NONSYMMETRIC  CODES, 

C  RESEARCH  REPORT  NO.  114,  DEPT.  OF  COMPUTER  SCIENCES,  YALE 

C  UNIVERSITY,  1977. 

C - 

C  SUMMARY  OF  USAGE. 

C 

C  COMMUNICATION  BETWEEN  THE  USER  AND  THE  LSODES  PACKAGE,  FOR  NORMAL 
C  SITUATIONS,  IS  SUMMARIZED  HERE.  THIS  SUMMARY  DESCRIBES  ONLY  A  SUBSET 
C  OF  THE  FULL  SET  OF  OPTIONS  AVAILABLE.  SEE  THE  FULL  DESCRIPTION  FOR 
C  DETAILS,  INCLUDING  OPTIONAL  COMMUNICATION,  NONSTANDARD  OPTIONS, 

C  AND  INSTRUCTIONS  FOR  SPECIAL  SITUATIONS.  SEE  ALSO  THE  EXAMPLE 
C  PROBLEM  (WITH  PROGRAM  AND  OUTPUT)  FOLLOWING  THIS  SUMMARY. 

C 

C  A.  FIRST  PROVIDE  A  SUBROUTINE  OF  THE  FORM.. 

C  SUBROUTINE  F  (NEQ,  T,  Y,  YDOT) 

C  DIMENSION  Y( NEQ) ,  YDOT (NEQ) 

C  WHICH  SUPPLIES  THE  VECTOR  FUNCTION  F  BY  LOADING  YDOT (I)  WITH  F(I) . 

C 

C  B.  NEXT  DETERMINE  (OR  GUESS)  WHETHER  OR  NOT  THE  PROBLEM  IS  STIFF. 


'  STIFFNESS  OCCURS  WHEN  THE  JACOBIAN  MATRIX  DF/DY  HAS  AN  EIGENVALUE 
C  WHOSE  REAL  PART  IS  NEGATIVE  AND  LARGE  IN  MAGNITUDE,  COMPARED  TO  THE 
^  RECIPROCAL  OF  THE  T  SPAN  OF  INTEREST.  IF  THE  PROBLEM  IS  NONSTIFF, 

USE  A  METHOD  FLAG  MF  »  10 .  IF  IT  IS  STIFF,  THERE  ARE  TWO-  STANDARD 
C  FOR  THE  METHOD  FLAG,  MF  =  121  AND  MF  =  222.  IN  BOTH  CASES,  LSODES 
_C  REQUIRES  THE  JACOBIAN  MATRIX  IN  SOME  FORM,  AND  IT  TREATS  THIS  MATRIX 
IN  GENERAL  SPARSE  FORM,  WITH  SPARSITY  STRUCTURE  DETERMINED  INTERNALLY 
(FOR  OPTIONS  WHERE  THE  USER  SUPPLIES  THE  SPARSITY  STRUCTURE,  SEE 
C  THE  FULL  DESCRIPTION  OF  MF  BELOW.) 

C.  IF  THE  PROBLEM  IS  STIFF,  YOU  ARE  ENCOURAGED  TO  SUPPLY  THE  JACOBIAN 
C  DIRECTLY  (MF  =  121) ,  BUT  IF  THIS  IS  NOT  FET^IBLE,  LSODES  WILL 
COMPUTE  IT  INTERNALLY  BY  DIFFERENCE  QUOTIENTS  (MF  =  222) . 

IF  YOU  ARE  SUPPLYING  THE  JACOBIAN,  PROVIDE  A  SUBROUTINE  OF  THE  FORM. . 
SUBROUTINE  JAC  (NEQ,  T,  Y,  J,  IAN,  JAN,  PDJ) 

C  DIMENSION  Y(l) ,  IAN (1) ,  JAN(l),  PDJ(l) 

HERE  NEQ,  T,  Y,  AND  J  ARE  INPUT  ARGUMENTS,  AND  THE  JAC  ROUTINE  IS  TO 
;  LOAD  THE  ARRAY  PDJ  (OF  LENGTH  NEQ)  WITH  THE  J-TH  COLUMN  OF  DF/DY. 

C  I.E.,  LOAD  PDJ (I)  WITH  DF(I)/DY(J)  FOR  ALL  RELEVANT  VALUES  OF  I. 

A  THE  ARGUMENTS  IAN  AND  JAN  SHOULD  BE  IGNORED  FOR  NORMAL  SITUATIONS. 

LSODES  WILL  CALL  THE  JAC  ROUTINE  WITH  J  =  1,2, .. . ,NEQ. 

V-  ONLY  NONZERO  ELEMENTS  NEED  BE  LOADED.  USUALLY,  A  CRUDE  APPROXIMATION 
C  TO  DF/DY,  POSSIBLY  WITH  FEWER  NONZERO . ELEMENTS ,  WILL  SUFFICE. 


D.  WRITE  A  MAIN  PROGRAM  WHICH  CALLS  SUBROUTINE  LSODES  ONCE  FOR 
C  EACH  POINT  AT  WHICH  ANSWERS  ARE  DESIRED.  THIS  SHOULD  ALSO  PROVIDE 
FOR  POSSIBLE  USE  OF  LOGICAL  UNIT  6  FOR  OUTPUT  OF  ERROR  MESSAGES 
BY  LSODES.  ON  THE  FIRST  CALL  TO  LSODES,  SUPPLY  ARGUMENTS  AS  FOLLOWS. 
C  F  =  NAME  OF  SUBROUTINE  FOR  RIGHT-HAND  SIDE  VECTOR  F. 

Q.  THIS  NAME  MUST  BE  DECLARED  EXTERNAL  IN  CALLING  PROGRAM. 

NEQ  =  NUMBER  OF  FIRST  ORDER  ODE-S. 

..  Y  =  ARRAY  OF  INITIAL  VALUES,  OF  LENGTH  NEQ. 

C  T  =  THE  INITIAL  VALUE  OF  THE  INDEPENDENT  VARIABLE. 

“  TOUT  =  FIRST  POINT  WHERE  OUTPUT  IS  DESIRED  (.NE.  T) . 

ITOL  =  1  OR  2  ACCORDING  AS  ATOL  (BELOW)  IS  A  SCALAR  OR  ARRAY. 

C  RTOL  =  RELATIVE  TOLERANCE  PARAMETER  (SCALAR) . 

^  ATOL  =  ABSOLUTE  TOLERANCE  PARAMETER  (SCALAR  OR  ARRAY) . 

THE  ESTIMATED  LOCAL  ERROR  IN  Y(I)  WILL  BE  CONTROLLED  SO  AS 
^  TO  BE  ROUGHLY  LESS  (IN  MAGNITUDE)  THAN 

C  EWT(I)  =  RTOL*ABS (Y(I) )  +  ATOL  IF  ITOL  =  1 ,  OR 

EWT(I)  =  RTOL*ABS (Y(I) )  +  ATOL (I)  IF  ITOL  =  2. 

THUS  THE  LOCAL  ERROR  TEST  PASSES  IF,  IN  EACH  COMPONENT, 

C  EITHER  THE  ABSOLUTE  ERROR  IS  LESS  THAN  ATOL  (OR  ATOL (I)) , 

OR  THE  RELATIVE  ERROR  IS  LESS  THAN  RTOL, 

USE  RTOL  =  0.0  FOR  PURE  ABSOLUTE  ERROR  CONTROL,  AND 
C  USE  ATOL  =0.0  (OR  ATOL ( I )  =  0.0)  FOR  PURE  RELATIVE  ERROR 

CL  CONTROL.  CAUTION. .  ACTUAL  (GLOBAL)  ERRORS  MAY  EXCEED  THESE 

LOCAL  TOLERANCES,  SO  CHOOSE  THEM  CONSERVATIVELY. 

ITASK  =  1  FOR  NORMAL  COMPUTATION  OF  OUTPUT  VALUES  OF  Y  AT  T  =  TOUT. 

C  ISTATE  =  INTEGER  FLAG  (INPUT  AND  OUTPUT) .  SET  ISTATE  =  1. 

lOPT  =  0  TO  INDICATE  NO  OPTIONAL  INPUTS  USED. 

RWORK  =  REAL  WORK  ARRAY  OF  LENGTH  AT  LEAST . . 

C  20  +  16*NEQ  FOR  MF  =  10, 

<5-  20  +  (2  +  l./LENRAT)  *NNZ  +  (11  +  9 . /LENRAT)  *NEQ 

FOR  MF  =  121  OR  222, 

WHERE .  . 

C,  NNZ  =  THE  NUMBER  OF  NONZERO  ELEMENTS  IN  THE  SPARSE 

JACOBIAN  (IF  THIS  IS  UNKNOWN,  USE  AN  ESTIMATE),  AND 
LENRAT  =  THE  REAL  TO  INTEGER  WORDLENGTH  RATIO  (USUALLY  1  IN 
C  SINGLE  PRECISION  AND  2  IN  DOUBLE  PRECISION) . 
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LRW 

IWORK 

LIW 

JAC 


MF 


INTEGER  WORK  ARRAY 
DECLARED  LENGTH  OF 
NAME  OF  SUBROUTINE 
IF  USED,  THIS  NAME 


IN  ANY  CASE,  THE  REQUIRED  SIZE  OF  RWORK  CANNOT  GENERALLY 
BE  PREDICTED  IN  ADVANCE  IF  MF  =  121  OR  222,  AND  THE  VALUE 
ABOVE  IS  A  ROUGH  ESTIMATE  OF  A  CRUDE  LOWER  BOUND.  SOME 
EXPERIMENTATION  WITH  THIS  SIZE  MAY  BE  NECESSARY. 

(WHEN  KNOWN,  THE  CORRECT  REQUIRED  LENGTH  IS  AN  OPTIONAL 
OUTPUT,  AVAILABLE  IN  IWORK (17).) 

DECLARED  LENGTH  OF  RWORK  (IN  USER-S  DIMENSION) . 

OF  LENGTH  AT  LEAST  30. 

IWORK  (IN  USER-S  DIMENSION) . 

FOR  UACOBIAN  MATRIX  (MF  =  121) . 

MUST  BE  DECIARED  EXTERNAL  IN  CALLING 
PROGRAM.  IF  NOT  USED,  PASS  A  DUMMY  NAME. 

METHOD  FLAG.  STANDARD  VALUES  ARE.. 

10  FOR  NONSTIFF  (ADAMS)  METHOD,  NO  JACOBIAN  USED. 

121  FOR  STIFF  (BDF)  METHOD,  USER-SUPPLIED  SPARSE  JACOBIAN. 
222  FOR  STIFF  METHOD,  INTERNALLY  GENERATED  SPARSE  JACOBIAN. 
NOTE  THAT  THE  MAIN  PROGRAM  MUST  DECLARE  ARRAYS  Y,  RWORK,  IWORK, 

AND  POSSIBLY  ATOL . 

E.  THE  OUTPUT  FROM  THE  FIRST  CALL  (OR  ANY  CALL)  IS.. 

Y  =  ARRAY  OF  COMPUTED  VALUES  OF  Y(T)  VECTOR. 

T  =  CORRESPONDING  VALUE  OF  INDEPENDENT  VARIABLE  (NORMALLY  TOUT) . 
ISTATE  =2  IF  LSODES  WAS  SUCCESSFUL,  NEGATIVE  OTHERWISE. 

-1  MEANS  EXCESS  WORK  DONE  ON  THIS  CALL  (PERHAPS  WRONG  MF) . 

-2  MEANS  EXCESS  ACCURACY  REQUESTED  (TOLERANCES  TOO  SMALL) . 

-3  MEANS  ILLEGAL  INPUT  DETECTED  (SEE  PRINTED  MESSAGE) . 

-4  MEANS  REPEATED  ERROR  TEST  FAILURES  (CHECK  ALL  INPUTS) . 

-5  MEANS  REPEATED  CONVERGENCE  FAILURES  (PERHAPS  BAD  JACOBIAN 
SUPPLIED  OR  WRONG  CHOICE  OF  MF  OR  TOLERANCES) . 

-6  MEANS  ERROR  WEIGHT  BECAME  ZERO  DURING  PROBLEM.  (SOLUTION 
COMPONENT  I  VANISHED,  AND  ATOL  OR  ATOL (I)  =0.) 

A  RETURN  WITH  ISTATE  =  -1,  -4,  OR  -5  MAY  RESULT  FROM  USING 
AN  INAPPROPRIATE  SPARSITY  STRUCTURE,  ONE  THAT  IS  QUITE 
DIFFERENT  FROM  THE  INITIAL  STRUCTURE.  CONSIDER  CALLING 
LSODES  AGAIN  WITH  ISTATE  =  3  TO  FORCE  THE  STRUCTURE  TO  BE 
REEVALUATED.  SEE  THE  FULL  DESCRIPTION  OF  ISTATE  BELOW. 

F.  TO  CONTINUE  THE  INTEGRATION  AFTER  A  SUCCESSFUL  RETURN,  SIMPLY 
RESET  TOUT  AND  CALL  LSODES  AGAIN .  NO  OTHER  PARAMETERS  NEED  BE  RESET . 


C  FULL  DESCRIPTION  OF  USER  INTERFACE  TO  LSODES. 

C 

C  THE  USER  INTERFACE  TO  LSODES  CONSISTS  OF  THE  FOLLOWING  PARTS. 

C 

C  I.  THE  CALL  SEQUENCE  TO  SUBROUTINE  LSODES,  WHICH  IS  A  DRIVER 
C  ROUTINE  FOR  THE  SOLVER.  THIS  INCLUDES  DESCRIPTIONS  OF  BOTH 

C  THE  CALL  SEQUENCE  ARGUMENTS  AND  OF  USER-SUPPLIED  ROUTINES. 

C  FOLLOWING  THESE  DESCRIPTIONS  IS  A  DESCRIPTION  OF 

C  OPTIONAL  INPUTS  AVAILABLE  THROUGH  THE  CALL  SEQUENCE,  AND  THEN 

C  A  DESCRIPTION  OF  OPTIONAL  OUTPUTS  (IN  THE  WORK  ARRAYS) . 

C 

C  II.  DESCRIPTIONS  OF  OTHER  ROUTINES  IN  THE  LSODES  PACKAGE  THAT  MAY  BE 
C  (OPTIONALLY)  CALLED  BY  THE  USER.  THESE  PROVIDE  THE  ABILITY  TO 

C  ALTER  ERROR  MESSAGE  HANDLING,  SAVE  AND  RESTORE  THE  INTERNAL 

C  COMMON,  AND  OBTAIN  SPECIFIED  DERIVATIVES  OF  THE  SOLUTION  Y(T) . 

C 

C  III.  DESCRIPTIONS  OF  COMMON  BLOCKS  TO  BE  DECLARED  IN  OVERLAY 
C  OR  SIMILAR  ENVIRONMENTS,  OR  TO  BE  SAVED  WHEN  DOING  AN  INTERRUPT 

C  OF  THE  PROBLEM  AND  CONTINUED  SOLUTION  LATER. 
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C  IV. 


DESCRIPTION  OF  TWO  SUBROUTINES  IN  THE  LSODES  PACKAGE,  EITHER  OF 
WHICH  THE  USER  MAY  REPLACE  WITH  HIS  OWN  VERSION,  IF  DESIRED. 
THESE  RELATE  TO  THE  MEASUREM©IT  OF  ERRORS. 


PART  I.  CALL  SEQUENCE. 

THE  CALL  SEQUENCE  PARAMETERS  USED  FOR  INPUT  ONLY  ARE 

F,  NEQ,  TOUT,  ITOL,  RTOL,  ATOL,  ITASK,  lOPT,  LRW,  LIW,  JAC,  MF, 
AND  THOSE  USED  FOR  BOTH  INPUT  AND  OUTPUT  ARE 
Y,  T,  ISTATE. 

THE  WORK  ARRAYS  RWORK  AND  IWORK  ARE  ALSO  USED  FOR  CONDITIONAL  AND 
OPTIONAL  INPUTS  AND  OPTIONAL  OUTPUTS.  (THE  TERM  OUTPUT  HERE  REFERS 
TO  THE  RETURN  FROM  SUBROUTINE  LSODES  TO  THE  USER-S  CALLING  PROGRAM.) 

THE  LEGALITY  OF  INPUT  PARAMETERS  WILL  BE  THOROUGHLY  CHECKED  ON  THE 
INITIAL  CALL  FOR  THE  PROBLEM,  BUT  NOT  CHECKED  THEREAFTER  UNLESS  A 
CHANGE  IN  INPUT  PARAMETERS  IS  FLAGGED  BY  ISTATE  =  3  ON  INPUT. 

THE  DESCRIPTIONS  OF  THE  CALL  ARGUMENTS  ARE  AS  FOLLOWS. 

F  =  THE  NAME  OF  THE  USER- SUPPLIED  SUBROUTINE  DEFINING  THE 

ODE  SYSTEM.  THE  SYSTEM  MUST  BE  PUT  IN  THE  FIRST-ORDER 
FORM  DY/DT  =  F(T,Y),  WHERE  F  IS  A  VECTOR-VALUED  FUNCTION 
OF  THE  SCALAR  T  AND  THE  VECTOR  Y.  SUBROUTINE  F  IS  TO 
COMPUTE  THE  FUNCTION  F.  IT  IS  TO  HAVE  THE  FORM 
SUBROUTINE  F  (NEQ,  T,  Y,  YDOT) 

DIMENSION  Y(l) ,  YDOT(l) 

WHERE  NEQ,  T,  AND  Y  ARE  INPUT,  AND  THE  ARRAY  YDOT  =  F(T,Y) 

IS  OUTPUT.  Y  AND  YDOT  ARE  ARRAYS  OF  LENGTH  NEQ. 

(IN  THE  DIMENSION  STATEMENT  ABOVE,  1  IS  A  DUMMY 
DIMENSION..  IT  CAN  BE  REPLACED  BY  ANY  VALUE.) 

SUBROUTINE  F  SHOULD  NOT  ALTER  Y(l) , . . . ,Y(NEQ) . 

F  MUST  BE  DECLARED  EXTERNAL  IN  THE  CALLING  PROGRAM. 

SUBROUTINE  F  MAY  ACCESS  USER-DEFINED  QUANTITIES  IN 
NEQ (2) , . . .  AND  Y(NEQ(1) +1) , . . .  IF  NEQ  IS  AN  ARRAY 
(DIMENSIONED  IN  F)  AND  Y  HAS  LENGTH  EXCEEDING  NEQ(l) . 

SEE  THE  DESCRIPTIONS  OF  NEQ  AND  Y  BELOW. 

NEQ  =  THE  SIZE  OF  THE  ODE  SYSTEM  (NUMBER  OF  FIRST  ORDER 

ORDINARY  DIFFERENTIAL  EQUATIONS) .  USED  ONLY  FOR  INPUT. 

NEQ  MAY  BE  DECREASED,  BUT  NOT  INCREASED,  DURING  THE  PROBLEM. 
IF  NEQ  IS  DECREASED  (WITH  ISTATE  =  3  ON  INPUT) ,  THE 
REMAINING  COMPONENTS  OF  Y  SHOULD  BE  LEFT  UNDISTURBED,  IF 
THESE  ARE  TO  BE  ACCESSED  IN  F  AND/OR  JAC. 

NORMALLY,  NEQ  IS  A  SCALAR,  AND  IT  IS  GENERALLY  REFERRED  TO 
AS  A  SCALAR  IN  THIS  USER  INTERFACE  DESCRIPTION.  HOWEVER, 

NEQ  MAY  BE  AN  ARRAY,  WITH  NEQ(l)  SET  TO  THE  SYSTEM  SIZE. 

(THE  LSODES  PACKAGE  ACCESSES  ONLY  NEQ(l) .)  IN  EITHER  CASE, 

C  THIS  PARAMETER  IS  PASSED  AS  THE  NEQ  ARGUMENT  IN  ALL  CALLS 

^  TO  F  AND  JAC.  HENCE,  IF  IT  IS  AN  ARRAY,  LOCATIONS 

NEQ (2) , . . .  MAY  BE  USED  TO  STORE  OTHER  INTEGER  DATA  AND'  PASS 
C"'  IT  TO  F  AND/OR  JAC.  SUBROUTINES  F  AND/OR  JAC  MUST  INCLUDE 

Q.  NEQ  IN  A  DIMENSION  STATEMENT  IN  THAT  CASE. 

L  Y  =  A  REAL  ARRAY  FOR  THE  VECTOR  OF  DEPENDENT  VARIABLES,  OF 

C  LENGTH  NEQ  OR  MORE.  USED  FOR  BOTH  INPUT  AND  OUTPUT  ON  THE 
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FIRST  CALL  (I STATE  =1),  AND  ONLY  FOR  OUTPUT  ON  OTHER  CALLS. 
•  ON  THE  FIRST  CALL,  Y  MUST  CONTAIN  THE  VECTOR  OF  INITIAL 
VALUES.  ON  OUTPUT,  Y  CONTAINS  THE  COMPUTED  SOLUTION  VECTOR, 
EVALUATED  AT  T.  IF  DESIRED,  THE  Y  ARRAY  MAY  BE  USED 
FOR  OTHER  PURPOSES  BETWEEN  CALLS  TO  THE  SOLVER. 

THIS  ARRAY  IS  PASSED  AS  THE  Y  ARGUMENT  IN  ALL  CALLS  TO 
F  AND  JAC.  HENCE  ITS  LENGTH  MAY  EXCEED  NEQ,  AND  LOCATIONS 
Y(NEQ+1) , . . .  MAY  BE  USED  TO  STORE  OTHER  REAL  DATA  AND 
PASS  IT  TO  F  AND/OR  JAC.  (THE  LSODES  PACKAGE  ACCESSES  ONLY 
Y(l)  ,  .  .  .  ,Y(NEQ)  .  ) 

=  THE  INDEPENDENT  VARIABLE.  ON  INPUT,  T  IS  USED  ONLY  ON  THE 
FIRST  CALL,  AS  THE  INITIAL  POINT  OF  THE  INTEGRATION. 

ON  OUTPUT,  AFTER- EACH  CALL,  T  IS  THE  VALUE  AT  WHICH  A 
COMPUTED  SOLUTION  Y  IS  EVALUATED  (USUALLY  THE  STU^E  AS  TOUT) . 
ON  AN  ERROR  RETURN,  T  IS  THE  FARTHEST  POINT  REACHED. 

=  THE  NEXT  VALUE  OF  T  AT  WHICH  A  COMPUTED  SOLUTION  IS  DESIRED. 
USED  ONLY  FOR  INPUT. 

WHEN  STARTING  THE  PROBLEM  (ISTATE  =1),  TOUT  MAY  BE  EQUAL 
TO  T  FOR  ONE  CALL,  THEN  SHOULD  .NE.  T  FOR  THE  NEXT  CALL. 

,  FOR  THE  INITIAL  T,  AN  INPUT  VALUE  OF  TOUT  .NE.  T  IS  USED 
IN  ORDER  TO  DETERMINE  THE  DIRECTION  OF  THE  INTEGRATION 
(I.E.  THE  ALGEBRAIC  SIGN  OF  THE  STEP  SIZES)  AND  THE  ROUGH 
SCALE  OF  THE  PROBLEM.  INTEGRATION  IN  EITHER  DIRECTION 
(FORWARD  OR  BACKWARD  IN  T)  IS  PERMITTED. 

IF  ITASK  =  2  OR  5  (ONE-STEP  MODES) ,  TOUT  IS  IGNORED  AFTER 
THE  FIRST  CALL  (I.E.  THE  FIRST  CALL  WITH  TOUT  . NE .  T) . 
OTHERWISE,  TOUT  IS  REQUIRED  ON  EVERY  CALL. 

IF  ITASK  =  1,  3,  OR  4,  THE  VALUES  OF  TOUT  NEED  NOT  BE 
MONOTONE,  BUT  A  VALUE  OF  TOUT  WHICH  BACKS  UP  IS  LIMITED 
TO  THE  CURRENT  INTERNAL  T  INTERVAL,  WHOSE  ENDPOINTS  ARE 
TCUR  -  HU  AND  TCUR  (SEE  OPTIONAL  OUTPUTS,  BELOW,  FOR 
TCUR  AND  HU) . 

=  AN  INDICATOR  FOR  THE  TYPE  OF  ERROR  CONTROL.  SEE 
DESCRIPTION  BELOW  UNDER  ATOL.  USED  ONLY  FOR  INPUT. 

=  A  RELATIVE  ERROR  TOLERANCE  PARAMETER,  EITHER  A  SCALAR  OR 
AN  ARRAY  OF  LENGTH  NEQ.  SEE  DESCRIPTION  BELOW  UNDER  ATOL. 
INPUT  ONLY. 

=  AN  ABSOLUTE  ERROR  TOLERANCE  PARAMETER,  EITHER  A  SCALAR  OR 
AN  ARRAY  OF  LENGTH  NEQ.  INPUT  ONLY. 

THE  INPUT  PARAMETERS  ITOL,  RTOL,  AND  ATOL  DETERMINE 
THE  ERROR  CONTROL  PERFORMED  BY  THE  SOLVER.  THE  SOLVER  WILL 
CONTROL  THE  VECTOR  E  =  (E(I))  OF  ESTIMATED  LOCAL  ERRORS 
IN  Y,  ACCORDING  TO  AN  INEQUALITY  OF  THE  FORM 

RMS-NORM  OF  (  E(I)/EWT(I)  )  . LE .  1, 

WHERE  EWT(I)  =  RTOL ( I ) *ABS ( Y ( I ) )  +  ATOL (I), 

AND  THE  RMS-NORM  (ROOT -MEAN -SQUARE  NORM)  HERE  IS 
RMS-NORM(V)  =  SQRT(SUM  V(I)**2  /  NEQ).  HERE  EWT  =  (EWT(I)) 
IS  A  VECTOR  OF  WEIGHTS  WHICH  MUST  ALWAYS  BE  POSITIVE,  AND 
THE  VALUES  OF  RTOL  AND  ATOL  SHOULD  ALL  BE  NON-NEGATIVE . 

THE  FOLLOWING  TABLE  GIVES  THE  TYPES  (SCALAR/ ARRAY)  OF 


c 


c 

S 

•i 


RTOL  AND  ATOL,  AND  THE  CORRESPONDING  FORM  OF  EWT(I) . 


ITOL  RTOL  ATOL 

1  SCALAR  SCALAR 

2  SCALAR  ARRAY 

3  ARRAY  SCALAR 

4  ARRAY  ARRAY 


EWT ( I ) 

RTOL*ABS ( Y ( I ) )  +  ATOL 
RTOL*ABS ( Y ( I ) )  +  ATOL ( I ) 
RTOL(I) *ABS(Y(I) )  +ATOL 
RTOL ( I ) *ABS ( Y ( I ) )  +  ATOL ( I ) 


C  WHEN  EITHER  OF  THESE  PARAMETERS  IS  A  SCALAR,  IT  NEED  NOT 

T  BE  DIMENSIONED  IN  THE  USER- S  CALLING  PROGRAM. 


C  IF  NONE  OF  THE  ABOVE  CHOICES  (WITH  ITOL,  RTOL,  AND  ATOL 

_C  FIXED  THROUGHOUT  THE  PROBLEM)  IS  SUITABLE,  MORE  GENERAL 

*  ERROR  CONTROLS:  CAN  BE  OBTAINED  BY  SUBSTITUTING 

USER-SUPPLIED  ROUTINES  FOR  THE  SETTING  OF  EWT  AND/OR  FOR 
C  THE  NORM  CALCULATION.  SEE  PART  IV  BELOW. 


IF  GLOBAL  ERRORS  ARE  TO  BE  ESTIMATED  BY  MAKING  A  REPEATED 
C  RUN  ON  THE  SAME  PROBLEM  WITH  SMALLER  TOLERANCES,  THEN  ALL 

-C  COMPONENTS  OF  RTOL  AND  ATOL  (I.E.  OF  EWT)  SHOULD  BE  SCALED 

DOWN  UNIFORMLY. 


C  ITASK 


C 


C 

r 

c 
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AN  INDEX  SPECIFYING  THE  TASK  TO  BE  PERFORMED. 

INPUT  ONLY.  ITASK  HAS  THE  FOLLOWING  VALUES  AND  MEANINGS. 

1  MEANS  NORMAL  COMPUTATION  OF  OUTPUT  VALUES  OF  Y{T)  AT 
T  =  TOUT  (BY  OVERSHOOTING  AND  INTERPOLATING) . 

2  MEANS  TAKE  ONE  STEP  ONLY  AND  RETURN. 

3  MEANS  STOP  AT  THE  FIRST  INTERNAL  MESH  POINT  AT  OR 
BEYOND  T  =  TOUT  AND  RETURN. 

4  MEANS  NORMAL  COMPUTATION  OF  OUTPUT  VALUES  OF  Y(T)  AT 
T  =  TOUT  BUT  WITHOUT  OVERSHOOTING  T  =  TCRIT. 

TCRIT  MUST  BE  INPUT  AS  RWORK(l) .  TCRIT  MAY  BE  EQUAL  TO 
OR  BEYOND  TOUT,  BUT  NOT  BEHIND  IT  IN  THE  DIRECTION  OF 
INTEGRATION.  THIS  OPTION  IS  USEFUL  IF  THE  PROBLEM 
HAS  A  SINGULARITY  AT  OR  BEYOND  T  =  TCRIT. 

5  MEANS  TAKE  ONE  STEP,  WITHOUT  PASSING  TCRIT,  AND  RETURN. 
TCRIT  MUST  BE  INPUT  AS  RWORK(l) . 


C 


C 

I  STATE 
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c 

A 
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NOTE. .  IF  ITASK  =  4  OR  5  AND  THE  SOLVER  REACHES  TCRIT 
(WITHIN  ROUNDOFF) ,  IT  WILL  RETURN  T  =  TCRIT  (EXACTLY)  TO 
INDICATE  THIS  (UNLESS  ITASK  =  4  AND  TOUT  COMES  BEFORE  TCRIT, 
IN  WHICH  CASE  ANSWERS  AT  T  =  TOUT  ARE  RETURNED  FIRST) . 

AN  INDEX  USED  FOR  INPUT  AND  OUTPUT  TO  SPECIFY  THE 
THE  STATE  OF  THE  CALCULATION. 

ON  INPUT,  THE  VALUES  OF  I STATE  ARE  AS  FOLLOWS. 

1  MEANS  THIS  IS  THE  FIRST  CALL  FOR  THE  PROBLEM 
(INITIALIZATIONS  WILL  BE  DONE) .  SEE  NOTE  BELOW. 

2  MEANS  THIS  IS  NOT  THE  FIRST  CALL,  AND  THE  CALCULATION 
IS  TO  CONTINUE  NORMALLY,  WITH  NO  CHANGE  IN  ANY  INPUT 
PARAMETERS  EXCEPT  POSSIBLY  TOUT  AND  ITASK. 

(IF  ITOL,  RTOL,  AND/OR  ATOL  ARE  CHANGED  BETWEEN  CALLS 
WITH  ISTATE  *  2,  THE  NEW  VALUES  WILL  BE  USED  BUT  NOT 
TESTED  FOR  LEGALITY.) 

3  MEANS  THIS  IS  NOT  THE  FIRST  CALL,  AND  THE 
CALCULATION  IS  TO  CONTINUE  NORMALLY,  BUT  WITH 
A  CHANGE  IN  INPUT  PARAMETERS  OTHER  THAN 

TOUT  AND  ITASK .  CHANGES  ARE  ALLOWED  IN 
NEQ,  ITOL,  RTOL,  ATOL,  lOPT,  LRW,  LIW,  MF, 
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THE  CONDITIONAL  INPUTS  lA  AND  JA, 

AND  ANY  OF  THE  OPTIONAL  INPUTS  EXCEPT  HO. 

IN  PARTICULAR,  IF  MITER  =  1  OR  2,  A  CALL  WITH  ISTATE  =  3 
WILL  CAUSE  THE  SPARSITY  STRUCTURE  OF  THE  PROBLEM  TO  BE 
RECOMPUTED  (OR  REREAD  FROM  I A  AND  JA  IF  MOSS  =  0) . 

NOTE . .  A  PRELIMINARY  CALL  WITH  TOUT  «  T  IS  NOT  COUNTED 
AS  A  FIRST  CALL  HERE,  AS  NO  INITIALIZATION  OR  CHECKING  OF 
INPUT  IS  DONE.  (SUCH  A  CiUjL  IS  SOMETIMES  USEFUL  FOR  THE 
PURPOSE  OF  OUTPUTTING  THE  INITIAL  CONDITIONS . ) 

THUS  THE  FIRST  CALL  FOR' WHICH  TOUT  .NE  .  T  REQUIRES 
ISTATE  =  1  ON  INPUT. 

ON  OUTPUT,  ISTATE  HAS  THE  FOLLOWING  VALUES  AND  MEANINGS. 

1  MEANS  NOTHING  WAS  DONE,  AS  TOUT  WAS  EQUAL  TO  T  WITH 
ISTATE  «  1  ON  INPUT.  (HOWEVER,  AN  INTERNAL  COUNTER  WAS 
SET  TO  DETECT  AND  PREVENT  REPEATED  CALLS  OF  THIS  TYPE.) 

2  MEANS  THE  INTEGRATION  WAS  PERFORMED  SUCCESSFULLY. 

-1  MEANS  AN  EXCESSIVE  AMOUNT  OF  WORK  (MORE  THAN  MXSTEP 

STEPS)  WAS  DONE  ON  THIS  CALL,  BEFORE  COMPLETING  THE 
REQUESTED  TASK,  BUT  THE  INTEGRATION  WAS  OTHERWISE 
SUCCESSFUL  AS  FAR  AS  T.  (MXSTEP  IS  AN  OPTIONAL  INPUT 
AND  IS  NORMALLY  500.)  TO  CONTINUE,  THE  USER  MAY 
SIMPLY  RESET  ISTATE  TO  A  VALUE  .GT.  1  AND  CALL  AGAIN 
(THE  EXCESS  WORK  STEP  COUNTER  WILL  BE  RESET  TO  0) . 

IN  ADDITION,  THE  USER  MAY  INCREASE  MXSTEP  TO  AVOID 
THIS  ERROR  RETURN  (SEE  BELOW  ON  OPTIONAL  INPUTS) . 

-2  MEANS  TOO  MUCH  ACCURACY  WAS  REQUESTED  FOR  THE  PRECISION 
OF  THE  MACHINE  BEING  USED.  THIS  WAS  DETECTED  BEFORE 
COMPLETING  THE  REQUESTED  TASK,  BUT  THE  INTEGRATION 
WAS  SUCCESSFUL  AS  FAR  AS  T.  TO  CONTINUE,  THE  TOLERANCE 
PARAMETERS  MUST  BE  RESET,  AND  ISTATE  MUST  BE  SET 
TO  3.  THE  OPTIONAL  OUTPUT  TOLSF  MAY  BE  USED  FOR  THIS 
PURPOSE.  (NOTE..  IF  THIS  CONDITION  IS  DETECTED  BEFORE 
TAKING  ANY  STEPS,  THEN  AN  ILLEGAL  INPUT  RETURN 
(ISTATE  =  -3)  OCCURS  INSTEAD.) 

-3  MEANS  ILLEGAL  INPUT  WAS  DETECTED,  BEFORE  TAKING  ANY 
INTEGRATION  STEPS.  SEE  WRITTEN  MESSAGE  FOR  DETAILS. 
NOTE..  IF  THE  SOLVER  DETECTS  AN  INFINITE  LOOP  OF  CALLS 
TO  THE  SOLVER  WITH  ILLEGAL  INPUT,  IT  WILL  CAUSE 
THE  RUN  TO  STOP. 

-4  MEANS  THERE  WERE  REPEATED  ERROR  TEST  FAILURES  ON 

ONE  ATTEMPTED  STEP,  BEFORE  COMPLETING  THE  REQUESTED 
TASK,  BUT  THE  INTEGRATION  WAS  SUCCESSFUL  AS  FAR  AS  T. 
THE  PROBLEM  MAY  HAVE  A  SINGULARITY,  OR  THE  INPUT 
MAY  BE  INAPPROPRIATE. 

-5  MEANS  THERE  WERE  REPEATED  CONVERGENCE  TEST  FAILURES  ON 
ONE  ATTEMPTED  STEP,  BEFORE  COMPLETING  THE  REQUESTED 
TASK,  BUT  THE  INTEGRATION  WAS  SUCCESSFUL  AS  FAR  AS  T. 
THIS  MAY  BE  CAUSED  BY  AN  INACCURATE  JACOBIAN  MATRIX, 

IF  ONE  IS.  BEING  USED. 

-6  MEANS  EWT(I)  BECAME  ZERO  FOR  SOME  I  DURING  THE 

INTEGRATION.  PURE  RELATIVE  ERROR  CONTROL  (ATOL ( I ) =0 . 0 ) 
WAS  REQUESTED  ON  A  VARIABLE  WHICH  HAS  NOW  VANISHED. 

THE  INTEGRATION  WAS  SUCCESSFUL  AS  FAR  AS  T . 

NOTE..  AN  ERROR  RETURN  WITH  ISTATE  =  -1,'  -4,  OR  -5  AND  WITH 
MITER  =  1  OR  2  MAY  MEAN  THAT  THE  SPARSITY  STRUCTURE  OF  THE 
PROBLEM  HAS  CHANGED  SIGNIFICANTLY  SINCE  IT  WAS  LAST 
DETERMINED  (OR  INPUT) .  IN  THAT  CASE,  ONE  CAN  ATTEMPT  TO 
COMPLETE  THE  INTEGRATION  BY  SETTING  ISTATE  =  3  ON  THE  NEXT 


CALL,  SO  THAT  A  NEW  STRUCTURE  DETERMINATION  IS  DONE. 


NOTE. .  SINCE  THE  NORMAL  OUTPUT  VALUE  OF  I STATE  IS  2, 

IT  DOES  NOT  NEED  TO  BE  RESET  FOR  NORMAL  CONTINUATION. 
ALSO,  SINCE  A  NEGATIVE  INPUT  VALUE  OF  ISTATE  WILL  BE 
REGARDED  AS  ILLEGAL,  A  NEGATIVE  OUTPUT  VALt®  REQUIRES  THE 
USER  TO  CHANGE  IT,  AND  POSSIBLY  OTHER  INPUTS,  BEFORE 
CALLING  THE  SOLVER  AGAIN. 


-T  lOPT  =  AN  INTEGER  FLAG  TO  SPECIFY  WHETHER  OR  NOT  ANY  OPTIONAL 
;  INPUTS  ARE  BEING  USED  ON  THIS  CALL.  INPUT  ONLY. 

C  THE  OPTIONAL  INPUTS  ARE  LISTED  SEPARATELY  BELOW. 

X  lOPT  =  0  MEANS  NO  OPTIONAL  INPUTS  ARE  BEING  USED. 

'  DEFAULT  VALUES  WILL  BE  USED  IN  ALL  CASES. 

lOPT  =  1  MEANS  ONE  OR  MORE  OPTIONAL  INPUTS  ARE  BEING  USED. 


C 

C 


RWORK  =  A  WORK  ARRAY  USED  FOR  A  MIXTURE  OF  REAL  (SINGLE  PRECISION) 
AND  INTEGER  WORK  SPACE. 

THE  LENGTH  OF  RWORK  (IN  REAL  WORDS)  MUST  BE  AT  LEAST 
20  +  NYH* (MAXORD  +  1)  +  3*NEQ  +  LWM  WHERE 
NYH  =  THE  INITIAL  VALUE  OF  NEQ, 

MAXORD  =12  (IF  METH  =  1)  OR  5  (IF  METH  =  2)  (UNLESS  A 
SMALLER  VALUE  IS  GIVEN  AS  AN  OPTIONAL  INPUT) , 


C 

X 


C 
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(NNZ+9*NEQ) /LENRAT 
(NNZ+10*NEQ) /LENRAT 


IP  MITER 
IF  MITER 
IF  MITER 
IF  MITER 


0, 
1, 
2, 
3  . 


LWM  =  0 

LWM  =  2*NNZ  +  2*NEQ  + 

LWM  =  2*NNZ  +  2*NEQ  + 

LWM  =  NEQ  +  2 
IN  THE  ABOVE  FORMULAS, 

NNZ  =  NUMBER  OF  NONZERO  ELEMENTS  IN  THE  JACOBIAN  MATRIX. 
LENRAT  =  THE  REAL  TO  INTEGER  WORDLENGTH  RATIO  (USUALLY  1  IN 
SINGLE  PRECISION  AND  2  IN  DOUBLE  PRECISION) . 

(SEE  THE  MF  DESCRIPTION  FOR  METH  AND  MITER.) 

THUS  IF  MAXORD  HAS  ITS  DEFAULT  VALUE  AND  NEQ  IS  CONSTANT, 
THE  MINIMUM  LENGTH  OF  RWORK  IS.. 


20 

+ 

16*NEQ 

FOR 

MF  =  10, 

c 

20 

+ 

16*NEQ  + 

LWM 

FOR 

MF  =11, 

Ill, 

211, 

12, 

112  , 

212, 

Xi. 

22 

+ 

17*NEQ 

FOR 

MF  =  13, 

20 

+ 

9*  NEQ 

FOR 

MF  =  20, 

20 

+ 

9*NEQ  + 

LWM 

FOR 

MF  =  21, 

121, 

221, 

22, 

122, 

222, 

_c 

22 

+ 

10*NEQ 

FOR 

MF  =  23. 

IF  MITER  =  1  OR  2,  THE  ABOVE  FORMULA  FOR  LWM  IS  ONLY  A 
CRUDE  LOWER  BOUND.  THE  REQUIRED  LENGTH  OF  RWORK  CANNOT 
BE  READILY  PREDICTED  IN  GENERAL,  AS  IT  DEPENDS  ON  THE 
SPARSITY  STRUCTURE  OF  THE  PROBLEM.  SOME  EXPERIMENTATION 
MAY  BE  NECESSARY. 

THE  FIRST  20  WORDS  OF  RWORK  ARE  RESERVED  FOR  CONDITIONAL 
AND  OPTIONAL  INPUTS  AND  OPTIONAL  OUTPUTS. 


C  THE  FOLLOWING  WORD  IN  RWORK  IS  A  CONDITIONAL  INPUT. . 

RWORK (1)  =  TCRIT  =  CRITICAL  VALUE  OF  T  WHICH  THE  SOLVER 
IS  NOT  TO  OVERSHOOT.  REQUIRED  IF  ITASK  IS 
C  4  OR  5,  AND  IGNORED  OTHERWISE.  (SEE  ITASK.) 

LRW  =  THE  LENGTH  OF  THE  ARRAY  RWORK,  AS  DECLARED  BY  THE  USER. 

C"  (THIS  WILL  BE  CHECKED  BY  THE  SOLVER.) 

C, 

IWORK  =  AN  INTEGER  WORK  ARRAY.  THE  LENGTH  OF  IWORK  MUST  BE  AT  LEAST 
31  +  NEQ  +  NNZ  IF  MOSS  =  0  AND  MITER  =  1  OR  2 ,  OR 
C  30  OTHERWISE. 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


LIW 


NOTE, 


(NNZ  IS  THE  NUMBER  OF  NONZERO  ELEMENTS  IN  DF/DY.) 

IN  LSODES,  IWORK  IS  USED  ONLY  FOR  CONDITIONAL  AND 
OPTIONAL  INPUTS  AND  OPTIONAL  OUTPUTS. 

THE  FOLLOWING  TWO  BLOCKS  OF  WORDS  IN  IWORK  ARE  CONDITIONAL 
INPUTS,  REQUIRED  IF  MOSS  =  0  AND  MITER  =  1  OR  2 ,  BUT  NOT 
OTHERWISE  (SEE  THE  DESCRIPTION  OF  MF  FOR  MOSS) . 

IWORK(30+J)  =  IA(J)  (J=l, . . . ,NEQ+1) 

IWORK(31+NEQ+K)  =  JACK)  (K=l ,  .  .  .  ,  NNZ) 

THE  TWO  ARRAYS  lA  AND  JA  DESCRIBE  THE  SPARSITY  STRUCTURE 
TO  BE  ASSUMED  FOR  THE  JACOBIAN  MATRIX.  JA  CONTAINS  THE  ROW 
INDICES  WHERE  NONZERO  ELEMENTS  OCCUR,  READING  IN  COLUMNWISE 
ORDER,  AND  lA  CONTAINS  THE  STARTING  LOCATIONS  IN  JA  OF  THE 
DESCRIPTIONS  OF  COLUMNS  1,...,NEQ,  IN  THAT  ORDER,  WITH 
lA ( 1 )  =  1 .  THUS ,  FOR  EACH  COLUMN  INDEX  J  =  1 , . . . , NEQ ,  THE 
VALUES  OF  THE  ROW  INDEX  I  IN  COLUMN  J  WHERE  A  NONZERO 
ELEMENT  MAY  OCCUR  ARE  GIVEN  BY 

I  =  JA(K) ,  WHERE  IA(J)  .LE.  K  .LT.  IA(J+1). 

IF  NNZ  IS  THE  TOTAL  NUMBER  OF  NONZERO  LOCATIONS  ASSUMED, 
THEN  THE  LENGTH  OF  THE  JA  ARRAY  IS  NNZ,  AND  IA{NEQ+1)  MUST 
BE  NNZ  +  1.  DUPLICATE  ENTRIES  ARE  NOT  ALLOWED. 

THE  LENGTH  OF  THE  ARRAY  IWORK,  AS  DECLARED  BY  THE  USER. 
(THIS  WILL  BE  CHECKED  BY  THE  SOLVER.) 


THE  WORK  ARRAYS  MUST  NOT  BE  ALTERED  BETWEEN  CALLS  TO  LSODES 
FOR  THE  SAME  PROBLEM,  EXCEPT  POSSIBLY  FOR  THE  CONDITIONAL  AND 
OPTIONAL  INPUTS,  AND  EXCEPT  FOR  THE  LAST  3*NEQ  WORDS  OF  RWORK. 

THE  LATTER  SPACE  IS  USED  FOR  INTERNAL  SCRATCH  SPACE,  AND  SO  IS 
AVAILABLE  FOR  USE  BY  THE  USER  OUTSIDE  LSODES  BETWEEN  CALLS,  IF 
DESIRED  (BUT  NOT  FOR  USE  BY  F  OR  JAC) . 

JAC  =  NAME  OF  USER-SUPPLIED  ROUTINE  (MITER  =  1  OR  MOSS  =  1)  TO 
COMPUTE  THE  JACOBIAN  MATRIX,  DF/DY,  AS  A  FUNCTION  OF 
THE  SCALAR  T  AND  THE  VECTOR  Y.  IT  IS  TO  HAVE  THE  FORM 
SUBROUTINE  JAC  (NEQ,  T,  Y,  J,  IAN,  JAN,  PDJ) 

DIMENSION  Y(l) ,  lAN(l),  JAN(l),  PDJ(l) 

WHERE  NEQ,  T,  Y,  J,  IAN,  AND  JAN  ARE  INPUT,  AND  THE  ARRAY 
PDJ,  OF  LENGTH  NEQ,  IS  TO  BE  LOADED  WITH  COLUMN  J 
OF  THE  JACOBIAN  ON  OUTPUT.  THUS  DF(I)/DY(J)  IS  TO  BE 
LOADED  INTO  PD J ( I )  FOR  ALL  RELEVANT  VALUES  OF  I. 

HERE  T  AND  Y  HAVE  THE  SAME  MEANING  AS  IN  SUBROUTINE  F, 

AND  J  IS  A  COLUMN  INDEX  (1  TO  NEQ) .  IAN  AND  JAN  ARE 
UNDEFINED  IN  CALLS  TO  JAC  FOR  STRUCTURE  DETERMINATION 
(MOSS  =  1) .  OTHERWISE,  IAN  AND  JAN  ARE  STRUCTURE 
DESCRIPTORS,  AS  DEFINED  UNDER  OPTIONAL  OUTPUTS  BELOW,  AND 
SO  CAN  BE  USED  TO  DETERMINE  THE  RELEVANT  ROW  INDICES  I,  IF 
DESIRED.  (IN  THE  DIMENSION  STATEMENT  ABOVE,  1  IS  A 
DUMMY  DIMENSION. .  IT  CAN  BE  REPLACED  BY  ANY  VALUE.) 

JAC  NEED  NOT  PROVIDE  DF/DY  EXACTLY.  A  CRUDE 
APPROXIMATION  (POSSIBLY  WITH  GREATER  SPARSITY)  WILL  DO. 

IN  ANY  CASE,  PDJ  IS  PRESET  TO  ZERO  BY  THE  SOLVER, 

SO  THAT  ONLY  THE  NONZERO  ELEMENTS  NEED  BE  LOADED  BY  JAC. 
CALLS  TO  JAC  ARE  MADE  WITH  J  =  1, . . . ,NEQ,  IN  THAT  ORDER,  AND 
EACH  SUCH  SET  OF  CALLS  IS  PRECEDED  BY  A  CALL  TO  F  WITH  THE 
SAME  ARGUMENTS  NEQ,  T,  AND  Y.  THUS  TO  GAIN  SOME  EFFICIENCY, 
INTERMEDIATE  QUANTITIES  SHARED  BY  BOTH  CALCULATIONS  MAY  BE 
SAVED  IN  A  USER  COMMON  BLOCK  BY  F  AND  NOT  RECOMPUTED  BY  JAC, 
IF  DESIRED.  JAC  MUST  NOT  ALTER  ITS  INPUT  ARGUMENTS. 
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JAC  MUST  BE  DECLARED  EXTERNAL  IN  THE  CALLING  PROGRAM. 

SUBROUTINE  JAC  MAY  ACCESS  USER-DEFINED  QUANTITIES  IN 
NEQ{2),...  AND  Y(NEQ{1) +1)  ,  .  .  .  IF.  NEQ  IS  MI  ARRAY 
(DIMENSIONED  IN  JAC)  AND  Y  HAS  LENGTH  EXCEEDING  NEQ(l) . 

SEE  THE  DESCRIPTIONS  OF  NEQ  AND  Y  ABOVE. 

THE  METHOD  FLAG.  USED  ONLY  FOR  INPUT. 

MF  HAS  THREE  DECIMAL  DIGITS--  MOSS,  METH,  MITER- - 
MF  =  100*MOSS  +  10*METH  +  MITER. 

MOSS  INDICATES  THE  METHOD  TO  BE  USED  TO  OBTAIN  THE  SPARSITY 
STRUCTURE  OF  THE  JACOBIAN  MATRIX  IF  MITER  =  1  OR  2 . . 

MOSS  =  0  MEANS  THE  USER  HAS  SUPPLIED  lA  AND  JA 
(SEE  DESCRIPTIONS  UNDER  I WORK  ABOVE) . 

MOSS  =  1  MEANS  THE  USER  HAS  SUPPLIED  JAC  (SEE  BELOW) 

AND  THE  STRUCTURE  WILL  BE  OBTAINED  FROM  NEQ 
INITIAL  CALLS  TO  JAC. 

MOSS  =  2  MEANS  THE  STOUCTURE  WILL  BE  OBTAINED  FROM  NEQ+1 
INITIAL  CALLS  TO  F. 

METH  INDICATES  THE  BASIC  LINEAR  MULTISTEP  METHOD. . 

METH  =  1  MEANS  THE  IMPLICIT  ADAMS  METHOD. 

METH  =  2  MEANS  THE  METHOD  BASED  ON  BACKWARD 
DIFFERENTIATION  FORMULAS  (BDF-S) . 

MITER  INDICATES  THE  CORRECTOR  ITERATION  METHOD. . 

MITER  =  0  MEANS  FUNCTIONAL  ITERATION  (NO  JACOBIAN  MATRIX 
IS  INVOLVED) . 

MITER  =  1  MEANS  CHORD  ITERATION  WITH  A  USER-SUPPLIED 
SPARSE  JACOBIAN,  GIVEN  BY  SUBROUTINE  JAC. 

MITER  =  2  MEANS  CHORD  ITERATION  WITH  AN  INTERNALLY 

GENERATED  (DIFFERENCE  QUOTIENT)  SPARSE  JACOBIAN 
(USING  NGP  EXTRA  CALLS  TO  F  PER  DF/DY  VALUE, 

WHERE  NGP  IS  AN  OPTIONAL  OUTPUT  DESCRIBED  BELOW . ) 
MITER  =  3  MEANS  CHORD  ITERATION  WITH  AN  INTERNALLY 

GENERATED  DIAGONAL  JACOBIAN  APPROXIMATION. 

(USING  1  EXTRA  CALL  TO  F  PER  DF/DY  EVALUATION) . 

IF  MITER  =  1  OR  MOSS  =  1,  THE  USER  MUST  SUPPLY  A  SUBROUTINE 
JAC  (THE  NAME  IS  ARBITRARY)  AS  DESCRIBED  ABOVE  UNDER  JAC. 
OTHERWISE,  A  DUMMY  ARGUMENT  CAN  BE  USED. 


THE  STANDARD  CHOICES  FOR  MF  ARE.. 

MF  =  10  FOR  A  NONSTIFF  PROBLEM, 

MF  =  21  OR  22  FOR  A  STIFF  PROBLEM  WITH  lA/JA  SUPPLIED 
(21  IF  JAC  IS  SUPPLIED,  22  IF  NOT) , 

MF  =  121  FOR  A  STIFF  PROBLEM  WITH  JAC  SUPPLIED, 

BUT  NOT  lA/JA, 

MF  =  222  FOR  A  STIFF  PROBLEM  WITH  NEITHER  lA/JA  NOR 
JAC  SUPPLIED. 

THE  SPARSENESS  STRUCTURE  CAN  BE  CHANGED  DURING  THE 
PROBLEM  BY  MAKING  A  CALL  TO  LSODES  WITH  ISTATE  =  3. 


OPTIONAL  INPUTS. 

THE  FOLLOWING  IS  A  LIST  OF  THE  OPTIONAL  INPUTS  PROVIDED  FOR  IN  THE 
CALL  SEQUENCE.  (SEE  ALSO  PART  II.)  FOR  EACH  SUCH  INPUT  VARIABLE, 
THIS  TABLE  LISTS  ITS  NAME  AS  USED  IN  THIS  DOCUMENTATION,  ITS 
LOCATION  IN  THE  CALL  SEQUENCE,  ITS  MEANING,  AND  THE  DEFAULT  VALUE. 
THE  USE  OF  ANY  OF  THESE  INPUTS  REQUIRES  lOPT  =  1,  AND  IN  THAT 
CASE  ALL  OF  THESE  INPUTS  ARE  EXAMINED.  A  VALUE  OF  ZERO  FOR  ANY 
OF  THESE  OPTIONAL  INPUTS  WILL  CAUSE  THE  DEFAULT  VALUE  TO  BE  USED. 
THUS  TO  USE  A  SUBSET  OF  THE  OPTIONAL  INPUTS,  SIMPLY  PRELOAD 
C  LOCATIONS  5  TO  10  IN  RWORK  AND  IWORK  TO  0 . 0  AND  0  RESPECTIVELY,  AND 
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THEN  SET  THOSE  OF  INTEREST  TO  NONZERO  VALUES. 


NAME 

HO 

HMAX 

HMIN 


SETH 


LOCATION 


RWORK{5) 


RWORK ( 6 ) 


RWORK ( 7 ) 


RWORK ( 8 ) 


MAXORD  IWORK ( 5 ) 


MEANING  AND  DEFAULT  VALUE 

THE  STEP  SIZE  TO  BE  ATTEMPTED  ON  THE  FIRST  STEP. 
THE  DEFAULT  VALUE  IS  DETERMINED  BY  THE  SOLVER. 

THE  MAXIMUM  ABSOLUTE  STEP  SIZE  ALLOWED. 

THE  DEFAULT  VALUE  IS  INFINITE. 

THE  MINIMUM  ABSOLUTE  STEP  SIZE  ALLOWED. 

THE  DEFAULT  VALUE  IS  0.  (THIS  LOWER  BOUND  IS  NOT 
ENFORCED  ON  THE  FINAL  STEP  BEFORE  REACHING  TCRIT 
WHEN  ITASK  =  4  OR  5 . ) 

THE  ELEMENT  THRESHHOLD  FOR  SPARSITY  DETERMINATION 
WHEN  MOSS  =  1  OR  2.  IF  THE  ABSOLUTE  VALUE  OF 
AN  ESTIMATED  JACOBIAN  ELEMENT  IS  .LE.  SETH,  IT 
WILL  BE  ASSUMED  TO  BE  ABSENT  IN  THE  STRUCTURE. 

THE  DEFAULT  VALUE  OF  SETH  IS  0. 

THE  MAXIMUM  ORDER  TO  BE  ALLOWED.  THE  DEFAULT 
VALUE  IS  12  IF  METH  =  1,  AND  5  IF  METH  =  2. 

IF  MAXORD  EXCEEDS  THE  DEFAULT  VALUE,  IT  WILL 
BE  REDUCED  TO  THE  DEFAULT  VALUE. 

IF  MAXORD  IS  CHANGED  DURING  THE  PROBLEM,  IT  MAY 
CAUSE  THE  CURRENT  ORDER  TO  BE  REDUCED. 

MAXIMUM  NUMBER  OF  (INTERNALLY  DEFINED)  STEPS 
ALLOWED  DURING  ONE  CALL  TO  THE  SOLVER. 

THE  DEFAULT  VALUE  IS  500. 

MAXIMUM  NUMBER  OF  MESSAGES  PRINTED  (PER  PROBLEM) 
WARNING  THAT  T  +  H  =  T  ON  A  STEP  (H  =  STEP  SIZE) . 
THIS  MUST  BE  POSITIVE  TO  RESULT  IN  A  NON-DEFAULT 
VALUE.  THE  DEFAULT  VALUE  IS  10. 

OPTIONAL  OUTPUTS. 

AS  OPTIONAL  ADDITIONAL  OUTPUT  FROM  LSODES,  THE  VARIABLES  LISTED 
BELOW  ARE  QUANTITIES  RELATED  TO  THE  PERFORMANCE  OF  LSODES 
WHICH  ARE  AVAILABLE  TO  THE  USER.  THESE  ARE  COMMUNICATED  BY  WAY  OF 
THE  WORK  ARRAYS,  BUT  ALSO  HAVE  INTERNAL  MNEMONIC  NAMES  AS  SHOWN. 
EXCEPT  WHERE  STATED  OTHERWISE,  ALL  OF  THESE  OUTPUTS  ARE  DEFINED 
ON  ANY  SUCCESSFUL  RETURN  FROM  LSODES,  AND  ON  ANY  RETURN  WITH 
ISTATE  =  -1,  -2,  -4,  -5,  OR  -6.  ON  AN  ILLEGAL  INPUT  RETURN 
(ISTATE  =  -3),  THEY  WILL  BE  UNCHANGED  FROM  THEIR  EXISTING  VALUES 
(IF  ANY),  EXCEPT  POSSIBLY  FOR  TOLSF,  LENRW,  AND  LENIW. 

ON  ANY  ERROR  RETURN,  OUTPUTS  RELEVANT  TO  THE  ERROR  WILL  BE  DEFINED, 
AS  NOTED  BELOW. 

NAME  LOCATION  MEANING 

HU  RWORK (11)  THE  STEP  SIZE  IN  T  LAST  USED  (SUCCESSFULLY) .  ^ 

HCUR  RWORK (12)  THE  STEP  SIZE  TO  BE  ATTEMPTED  ON  THE  NEXT  STEP. 

TCUR  RWORK (13)  THE  CURRENT  VALUE  OF  THE  INDEPENDENT  VARIABLE 

WHICH  THE  SOLVER  HAS  ACTUALLY  REACHED,  I.E.  THE 
CURRENT  INTERNAL  MESH  POINT  IN  T.  ON  OUTPUT,  TCUR 


MXSTEP  IWORK (6) 


MXHNIL  IWORK (7) 


jo  Tv,  Jin  !o(  ho.y^n  lor  ho  Inc  b  c  b  n  .  Jo  r  ...  h  o  .  Jo(-,.ho 


WILL  ALWAYS  BE  AT  LEAST  AS  FAR  AS  THE  ARGUMENT 
T,  BUT  MAY  BE  FARTHER  (IF  INTERPOLATION  WAS  DONE) . 

RW0RK(14)  A  TOLERANCE  SCALE  FACTOR,  GREATER  THAN  1.0, 

COMPUTED  WHEN  A  REQUEST  FOR  TOO  MUCH  ACCURACY  WAS 
DETECTED  (ISTATE  =  -3  IF  DETECTED  AT  THE  START  OF 
THE  PROBLEM,  ISTATE  =  -2  OTHERWISE) .  IF  ITOL  IS 
LEFT  UNALTERED  BUT  RTOL  AND  ATOL  ARE  UNIFORMLY 
SCALED  UP  BY  A  FACTOR  OF  TOLSF  FOR  THE  NEXT  CALL, 
THEN  THE  SOLVER  IS  DEEMED  LIKELY  TO  SUCCEED. 

(THE  USER  MAY  ALSO  IGNORE  TOLSF  AND  ALTER  THE 
TOLERANCE  PARAMETERS  IN  ANY  OTHER  WAY  APPROPRIATE . ) 

IWORK(ll)  THE  NUMBER  OF  STEPS  TAKEN  FOR  THE  PROBLEM  SO  FAR. 

NFE  IWORK{12)  THE  NUMBER  OF  F  EVALUATIONS  FOR  THE  PROBLEM  SO  FAR, 

EXCLUDING  THOSE  FOR  STRUCTURE  DETERMINATION 
(MOSS  =  2) . 

NJE  IWORK(13)  THE  NUMBER  OF  JACOBIAN  EVALUATIONS  FOR  THE  PROBLEM 

SO  FAR,  EXCLUDING  THOSE  FOR  STRUCTURE  DETERMINATION 
(MOSS  =  1) . 

NQU  IWORK(14)  THE  METHOD  ORDER  LAST  USED  (SUCCESSFULLY) . 

NQCUR  IWORKdS)  THE  ORDER  TO  BE  ATTEMPTED  ON  THE  NEXT  STEP. 

IMXER  IWORK(16)  THE  INDEX  OF  THE  COMPONENT  OF  LARGEST  MAGNITUDE  IN 

THE  WEIGHTED  LOCAL  ERROR  VECTOR  (  E(I)/EWT(I)  ), 

ON  AN  ERROR  RETURN  WITH  ISTATE  =  -4  OR  -5. 

LENRW  IWORK(17)  THE  LENGTH  OF  RWORK  ACTUALLY  REQUIRED. 

THIS  IS  DEFINED  ON  NORMAL  RETURNS  AND  ON  AN  ILLEGAL 
INPUT  RETURN  FOR  INSUFFICIENT  STORAGE. 

LENIW  IWORKdS)  THE  LENGTH  OF  IWORK  ACTUALLY  REQUIRED. 

THIS  IS  DEFINED  ON  NORMAL  RETURNS  AND  ON  AN  ILLEGAL 
INPUT  RETURN  FOR  INSUFFICIENT  STORAGE. 

NNZ  IWORK (19)  THE  NUMBER  OF  NONZERO  ELEMENTS  IN  THE  JACOBIAN 

MATRIX,  INCLUDING  THE  DIAGONAL  (MITER  =  1  OR  2 ) . 
(THIS  MAY  DIFFER  FROM  THAT  GIVEN  BY  IA(NEQ+1)-1 
IF  MOSS  =  0,  BECAUSE  OF  ADDED  DIAGONAL  ENTRIES.) 

NGP  IWORK (20)  THE  NUMBER  OF  GROUPS  OF  COLUMN  INDICES,  USED  IN 

DIFFERENCE  QUOTIENT  JACOBIAN  APROXIMATIONS  IF 
MITER  =  2.  THIS  IS  ALSO  THE  NUMBER  OF  EXTRA  F 
EVALUATIONS  NEEDED  FOR  EACH  JACOBIAN  EVALUATION. 

NLU  IWORK (21)  THE  NUMBER  OF  SPARSE  LU  DECOMPOSITIONS  FOR  THE 

PROBLEM  SO  FAR. 

LYH  IWORK (22)  THE  BASE  ADDRESS  IN  RWORK  OF  THE  HISTORY  ARRAY  YH, 

DESCRIBED  BELOW  IN  THIS  LIST. 

IPIAN  IWORK (23)  THE  BASE  ADDRESS  OF  THE  STRUCTURE  DESCRIPTOR  ARRAY 

IAN,  DESCRIBED  BELOW  IN  THIS  LIST. 

IPJAN  IWORK (24)  THE  BASE  ADDRESS  OF  THE  STRUCTURE  DESCRIPTOR  ARRAY 
C  JAN,  DESCRIBED  BELOW  IN  THIS  LIST. 


C 

:  TOLSF 


c 


NST 


NZL 


NZU 


.IWORK(25) 


IWORK(26) 


THE  NUMBER  OF  NONZERO  ELEMENTS  IN  THE  STRICT  LOWER 
TRIANGLE  OF  THE  LU  FACTORIZATION  USED  IN  THE  CHORD 
ITERATION  (MITER  «  1  OR  2) . 

THE  NUMBER  OF  NONZERO  ELEMENTS  IN  THE  STRICT  UPPER 
TRIANGLE  OF  THE  LU  FACTORIZATION  USED  IN  THE  CHORD 
ITERATION  (MITER  =  1  OR  2) . 

THE  TOTAL  NUMBER  OF  NONZEROS  IN  THE  FACTORIZATION 
IS  THEREFORE  NZL  +  NZU  +  NEQ. 
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C  THE  FOLLOWING  FOUR  ARRAYS  ARE  SEC»1ENTS  OF  THE  RWORK  ARRAY  WHICH 
C  MAY  ALSO  BE  OF  INTEREST  TO  THE  USER  AS  OPTIONAL  OUTPUTS . 

FOR  EACH  ARRAY,  THE  TABLE  BELOW  GIVES  ITS  INTERNAL  NAME, 

ITS  BASE  ADDRESS,  AND  ITS  DESCRIPTION. 

FOR  YH  AND  ACOR,  THE  BASE  ADDRESSES  ARE  IN  RWORK  (A  REAL  ARRAY) . 

THE  INTEGER  ARRAYS  IAN  AND  JAN  ARE  TO  BE  OBTAINED  BY  DECLARING  AN 
INTEGER  ARRAY  IWK  AND  IDENTIFYING  IWK(l)  WITH  RWORK(21),  USING  EITHER 
AN  EQUIVALENCE  STATEMENT  OR  A  SUBROUTINE  CALL.  THEN  THE  BASE 
ADDRESSES  IPIAN  (OF  IAN)  AND  IPJAN  (OF  JAN)  IN  IWK  ARE  TO  BE  OBTAINED 
AS  OPTIONAL  OUTPUTS  IWORK(23)  AND  IWORK(24),  RESPECTIVELY. 

THUS  lAN(l)  IS  IWK (IPIAN),  ETC. 


NAME 

IAN 

JAN 


BASE  ADDRESS 

IPIAN  (IN  IWK) 
IPJAN  (IN  IWK) 
(SEE  ABOVE) 


DESCRIPTION 


+  1 


YH 


LYH 

(OPTIONAL 

OUTPUT) 


ACOR 


LENRW-NEQ+1 


STRUCTURE  DESCRIPTOR  ARRAY  OF  SIZE  NEQ 
STRUCTURE  DESCRIPTOR  ARRAY  OF  SIZE  NNZ . 

IAN  AND  JAN  TOGETHER  DESCRIBE  THE  SPARSITY 
STRUCTURE  OF  THE  JACOBIAN  MATRIX,  AS  USED  BY 
LSODES  WHEN  MITER  »  1  OR  2 . 

JAN  CONTAINS  THE  ROW  INDICES  OF  THE  NONZERO 
LOCATIONS,  READING  IN  COLUMNWISE  ORDER,  AND 
IAN  CONTAINS  THE  STARTING  LOCATIONS  IN  JAN  OF 
THE  DESCRIPTIONS  OF  COLUMNS  1, . . . ,NEQ,  IN 
THAT  ORDER,  WITH  lAN(l)  =  1.  THUS  FOR  EACH 
J  =  1,...,NEQ,  THE  ROW  INDICES  I  OF  THE 
NONZERO  LOCATIONS  IN  COLUMN  J  ARE 
I  =  JAN(K),  lAN(J)  .LE.  K  . LT .  IAN(J+1). 

NOTE  THAT  IAN(NEQ+1)  =  NNZ  +  1. 

(IF  MOSS  =  0,  IAN/ JAN  MAY  DIFFER  FROM  THE 
INPUT  lA/JA  BECAUSE  OF  A  DIFFERENT  ORDERING 
IN  EACH  COLUMN,  AND  ADDED  DIAGONAL  ENTRIES.) 

THE  NORDSIECK  HISTORY  ARRAY,  OF  SIZE  NYH  BY 
(NQCUR  +  1) ,  WHERE  NYH  IS  THE  INITIAL  VALUE 
OF  NEQ.  FOR  J  =  0 , 1 ,..., NQCUR,  COLUMN  J+1 
OF  YH  CONTAINS  HOUR** J/FACTORIAL ( J)  TIMES 
THE  J-TH  DERIVATIVE  OF  THE  INTERPOLATING 
POLYNOMIAL  CURRENTLY  REPRESENTING  THE  SOLUTION, 
EVALUATED  AT  T  =  TCUR .  THE  BASE  ADDRESS  LYH 
IS  ANOTHER  OPTIONAL  OUTPUT,  LISTED  ABOVE. 

ARRAY  OF  SIZE  NEQ  USED  FOR  THE  ACCUMULATED 
CORRECTIONS  ON  EACH  STEP,  SCALED  ON  OUTPUT 
TO  REPRESENT  THE  ESTIMATED  LOCAL  ERROR  IN  Y  - 
ON  THE  LAST  STEP.  THIS  IS  THE  VECTOR  E  IN 
THE  DESCRIPTION  OF  THE  ERROR  CONTROL.  IT  IS 
DEFINED  ONLY  ON  A  SUCCESSFUL  RETURN  FROM 
LSODES . 
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C  PART  II,  OTHER  ROUTINES  CALLABLE. 

THE  FOLLOWING  ARE  OPTIONAL  CALLS  WHICH  THE  USER  MAY  MAKE  TO 
C  GAIN  ADDITIONAL  CAPABILITIES  IN  CONJUNCTION  WITH  LSODES . 

(THE  ROUTINES  XSETUN  AND  XSETF  ARE  DESIGNED  TO  CONFORM  TO  THE 
SLATEC  ERROR  HANDLING  PACKAGE . ) 

C  FORM  OF  CALL  FUNCTION 

-n  CALL  XSETUN (LUN)  SET  THE  LOGICAL  UNIT  NUMBER,  LUN,  FOR 

OUTPUT  OF  MESSAGES  FROM  LSODES,  IF 
C  THE  DEFAULT  IS  NOT  DESIRED. 

jC  the  default  value  of  LUN  IS  6. 

CALL  XSETF (MFLAG)  SET  A  FLAG  TO  CONTROL  THE  PRINTING  OF 

MESSAGES  BY  LSODES. 

MFLAG  =  0  MEANS  DO  NOT  PRINT.  (DANGER.. 
THIS  RISKS  LOSING  VALUABLE  INFORMATION.) 
MFLAG  =  1  MEANS  PRINT  (THE  DEFAULT) . 

EITHER  OF  THE  ABOVE  CALLS  MAY  BE  MADE  AT 
ANY  TIME  AND  WILL  TAKE  EFFECT  IMMEDIATELY. 

CALL  SVCMS  (RSAV,  ISAV)  STORE  IN  RSAV  AND  ISAV  THE  CONTENTS 

OF  THE  INTERNAL  COMMON  BLOCKS  USED  BY 
LSODES  (SEE  PART  III  BELOW) . 

RSAV  MUST  BE  A  REAL  ARRAY  OF  LENGTH  225 
^  OR  MORE,  AND  ISAV  MUST  BE  AN  INTEGER 

ARRAY  OF  LENGTH  75  OR  MORE. 

CALL  RSCMS  (RSAV,  ISAV)  RESTORE,  FROM  RSAV  AND  ISAV,  THE  CONTENTS 

OF  THE  INTERNAL  COMMON  BLOCKS  USED  BY 
LSODES .  PRESUMES  A  PRIOR  CALL  TO  SVCMS 
WITH  THE  SAME  ARGUMENTS. 

SVCMS  AND  RSCMS  ARE  USEFUL  IF 
INTERRUPTING  A  RUN  AND  RESTARTING 
LATER,  OR  ALTERNATING  BETWEEN  TWO  OR 
MORE  PROBLEMS  SOLVED  WITH  LSODES. 

CALL  INTDY(, , , , , )  PROVIDE  DERIVATIVES  OF  Y,  OF  VARIOUS 

(SEE  BELOW)  ORDERS,  AT  A  SPECIFIED  POINT  T,  IF 

DESIRED.  IT  MAY  BE  CALLED  ONLY  AFTER 
A  SUCCESSFUL  RETURN  FROM  LSODES. 

THE  DETAILED  INSTRUCTIONS  FOR  USING  INTDY  ARE  AS  FOLLOWS. 

THE  FORM  OF  THE  CALL  IS. . 

LYH  =  IWORK(22) 

CALL  INTDY  (T,  K,  RWORK(LYH),  NYH,  DKY,  I FLAG) 

THE  INPUT  PARAMETERS  ARE . . 

C  . 

T  =  VALUE  OF  INDEPENDENT  VARIABLE  WHERE  ANSWERS  ARE  DESIRED 

^  (NORMALLY  THE  SAME  AS  THE  T  LAST  RETURNED  BY  LSODES) . 

C  FOR  VALID  RESULTS,  T  MUST  LIE  BETWEEN  TCUR  -  HU  AND  TCUR . 

a  (SEE  OPTIONAL  OUTPUTS  FOR  TCUR  AND  HU.) 

K  =  INTEGER  ORDER  OF  THE  DERIVATIVE  DESIRED.  K  MUST  SATISFY 

0  .LE.  K  .LE.  NQCUR,  WHERE  NQCUR  IS  THE  CURRENT  ORDER 
C  (SEE  OPTIONAL  OUTPUTS) .  THE  CAPABILITY  CORRESPONDING 
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LYH 


NYH 


TO  K  =  0,  I.E.  COMPUTING  Y(T),  IS  ALREADY  PROVIDED 
BY  LSODES  DIRECTLY.  SINCE  NQCUR  .GE.  1,  THE  FIRST 
DERIVATIVE  DY/DT  IS  ALWAYS  AVAILABLE  WITH  INTDY. 

THE  BASE  ADDRESS  OF  THE  HISTORY  ARRAY  YH,  OBTAINED 
AS  AN  OPTIONAL  OUTPUT  AS  SHOWN  ABOVE. 

COLUMN  LENGTH  OF  YH,  EQUAL  TO  THE  INITIAL  VALUE  OF  NEQ. 


THE  OUTPUT  PARAMETERS  ARE . . 

DKY  =  A  REAL  ARRAY  OF  LENGTH  NEQ  CONTAINING  THE  COMPUTED  VALUE 

OF  THE  K-TH. DERIVATIVE  OF  Y{T) . 

I FLAG  =  INTEGER  FLAG,  RETURNED  AS  0  IF  K  AND  T  WERE  LEGAL, 

-1  IF  K  WAS  ILLEGAL,  AND  -2  IF  T  WAS  ILLEGAL. 

ON  AN  ERROR  RETURN,  A  MESSAGE  IS  ALSO  WRITTEN. 


C  PART  III. 
C 


COMMON  BLOCKS, 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


IF  LSODES  IS  TO  BE  USED  IN  AN  OVERLAY  SITUATION,  THE  USER 
MUST  DECLARE,  IN  THE  PRIMARY  OVERLAY,  THE  VARIABLES  IN. . 

(1)  THE  CALL  SEQUENCE  TO  LSODES, 

(2)  THE  THREE  INTERNAL  COMMON  BLOCKS 

/LSOOOl/  OF  LENGTH  258  {219  SINGLE  PRECISION  WORDS 

FOLLOWED  BY  39  INTEGER  WORDS) , 
/LSSOOl/  OF  LENGTH  40  (  6 

FOLLOWED  BY  34 


/EHOOOl/  OF  LENGTH  2  (INTEGER  WORDS) 


SINGLE  PRECISION  WORDS 
INTEGER  WORDS) , 


C  IF  LSODES  IS  USED  ON  A  SYSTEM  IN  WHICH  THE  CONTENTS  OF  INTERNAL 
C  COMMON  BLOCKS  ARE  NOT  PRESERVED  BETWEEN  CALLS,  THE  USER  SHOULD 
C  DECLARE  THE  ABOVE  THREE  COMMON  BLOCKS  IN  HIS  MAIN  PROGRAM  TO  INSURE 
C  THAT  THEIR  CONTENTS  ARE  PRESERVED. 

C 


C  IF  THE  SOLUTION  OF  A  GIVEN  PROBLEM  BY  LSODES  IS  TO  BE  INTERRUPTED 
C  AND  THEN  LATER  CONTINUED,  SUCH  AS  WHEN  RESTARTING  AN  INTERRUPTED  RUN 
C  OR  ALTERNATING  BETWEEN  TWO  OR  MORE  PROBLEMS,  THE  USER  SHOULD  SAVE, 

C  FOLLOWING  THE  RETURN  FROM  THE  LAST  LSODES  CALL  PRIOR  TO  THE 
C  INTERRUPTION,  THE  CONTENTS  OF  THE  CALL  SEQUENCE  VARIABLES  AND  THE 
C  INTERNAL  COMMON  BLOCKS,  AND  LATER  RESTORE  THESE  VALUES  BEFORE  THE 
C  NEXT  LSODES  CALL  FOR  THAT  PROBLEM.  TO  SAVE  AND  RESTORE  THE  COMMON 
C  BLOCKS,  USE  SUBROUTINES  SVCMS  AND  RSCMS  (SEE  PART  II  ABOVE) . 

C 

C  NOTE..  IN  THIS  VERSION  OF  LSODES,  THERE  ARE  TWO  DATA  STATEMENTS, 

C  IN  SUBROUTINES  LSODES  AND  XERRWV,  WHICH  LOAD  VARIABLES  INTO  THESE 
C  LABELED  COMMON  BLOCKS.  ON  SOME  SYSTEMS,  IT  MAY  BE  NECESSARY  TO 
C  MOVE  THESE  TO  A  SEPARATE  BLOCK  DATA  SUBPROGRAM. 

C 

C - 

C  PART  IV.  OPTIONALLY  REPLACEABLE  SOLVER  ROUTINES. 

C 

C  BELOW  ARE  DESCRIPTIONS  OF  TWO  ROUTINES  IN  THE  LSODES  PACKAGE  WHICH 
C  RELATE  TO  THE  MEASUREMENT  OF  ERRORS.  EITHER  ROUTINE  CAN  BE 
C  REPLACED  BY  A  USER-SUPPLIED  VERSION,  IF  DESIRED.  HOWEVER,  SINCE  SUCH 
C  A  REPLACEMENT  MAY  HAVE  A  MAJOR  IMPACT  ON  PERFORMANCE,  IT  SHOULD  BE 
C  DONE  ONLY  WHEN  ABSOLUTELY  NECESSARY,  AND  ONLY  WITH  GREAT  CAUTION. 

C  (NOTE. .  THE  MEANS  BY  WHICH  THE  PACKAGE  VERSION  OF  A  ROUTINE  IS  ‘ 

C  SUPERSEDED  BY  THE  USER-S  VERSION  MAY  BE  SYSTEM - DEPENDENT . ) 

C 

C  (A)  EWSET. 

C  THE  FOLLOWING  SUBROUTINE  IS  CALLED  JUST  BEFORE  EACH  INTERNAL 
C  INTEGRATION  STEP,  AND  SETS  THE  ARRAY  OF  ERROR  WEIGHTS,  EWT,  AS 
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:  DESCRIBED  UNDER  ITOL/RTOL/ATOL  ABOVE.. 

C  SUBROUTINE  EWSET  (NEQ,  ITOL,  RTOL,  ATOL,  YCUR,  EWT) 

WHERE  NEQ,  ITOL,  RTOL,  AND  ATOL  ARE  AS  IN  THE  LSODES  CALL  SEQUENCE, 

:  YCUR  CONTAINS  THE  CURRENT  DEPENDENT  VARIABLE  VECTOR,  AND 
C  EWT  IS  THE  ARRAY  OF  WEIGHTS  SET  BY  EWSET. 

s: 

!  IF  THE  USER  SUPPLIES  THIS  SUBROUTINE,  IT  MUST  RETURN  IN  EWT (I) 

J  (I  =  1, . , .,NEQ)  A  POSITIVE  QUANTITY  SUITABLE  FOR  COMPARING  ERRORS 
C  IN  Y(I)  TO.  THE  EWT  ARRAY  RETURNED  BY  EWSET  IS  PASSED  TO  THE 
T  VNORM  ROUTINE  (SEE  BELOW),  AND  ALSO  USED  BY  LSODES  IN  THE  COMPUTATION 
'  OF  THE  OPTIONAL  OUTPUT  IMXER,  THE  DIAGONAL  JACOBIAN  APPROXIMATION, 

C  AND  THE  INCREMENTS  FOR  DIFFERENCE  QUOTIENT  JACOB IANS . 

-C 

i  IN  THE  USER- SUPPLIED  VERSION  OF  EWSET,  IT  MAY  BE  DESIRABLE  TO  USE 
THE  CURRENT  VALUES  OF  DERIVATIVES  OF  Y.  DERIVATIVES  UP  TO  ORDER  NQ 
ARE  AVAILABLE  FROM  THE  HISTORY  ARRAY  YH,  DESCRIBED  ABOVE  UNDER 
OPTIONAL  OUTPUTS.  IN  EWSET,  YH  IS  IDENTICAL  TO  THE  YCUR  ARRAY, 
EXTENDED  TO  NQ  +  1  COLUMNS  WITH  A  COLUNBJ  LENGTH  OF  NYH  AND  SCALE 
FACTORS  OF  H* * J/FACTORIAL ( J) .  ON  THE  FIRST  CALL  FOR  THE  PROBLEM, 
GIVEN  BY  NST  =  0,  NQ  IS  1  AND  H  IS  TEMPORARILY  SET  TO  1.0. 

THE  QUANTITIES  NQ,  NYH,  H,  AND  NST  CAN  BE  OBTAINED  BY  INCLUDING 
IN  EWSET  THE  STATEMENTS . . 

COMMON  /LSOOOl/  RLS (219) , ILS (39) 

NQ  =  ILS (35) 

NYH  =  ILS (14) 

NST  =  ILS (36) 

H  =  RLS (213) 

THUS,  FOR  EXAMPLE,  THE  CURRENT  VALUE  OF  DY/DT  CAN  BE  OBTAINED  AS 
YCUR(NYH+I) /H  ( 1=1 , . . . , NEQ)  (AND  THE  DIVISION  BY  H  IS 
UNNECESSARY  WHEN  NST  =  0) . 

(B)  VNORM. 

THE  FOLLOWING  IS  A  REAL  FUNCTION  ROUTINE  WHICH  COMPUTES  THE  WEIGHTED 
ROOT -MEAN -SQUARE  NORM  OF  A  VECTOR  V. . 

D  =  VNORM  (N,  V,  W) 

WHERE . . 

N  =  THE  LENGTH  OF  THE  VECTOR, 

V  =  REAL  ARRAY  OF  LENGTH  N  CONTAINING  THE  VECTOR, 

W  =  REAL  ARRAY  OF  LENGTH  N  CONTAINING  WEIGHTS, 

D  =  SQRT(  (1/N)  *  SUM(V(I) *W(I) ) **2  ). 

VNORM  IS  CALLED  WITH  N  =  NEQ  AND  WITH  W(I)  =  1.0/EWT(I),  WHERE 
EWT  IS  AS  SET  BY  SUBROUTINE  EWSET. 

IF  THE  USER  SUPPLIES  THIS  FUNCTION,  IT  SHOULD  RETURN  A  NON-NEGATIVE 
VALUE  OF  VNORM  SUITABLE  FOR  USE  IN  THE  ERROR  CONTROL  IN  LSODES. 

NONE  OF  THE  ARGUMENTS  SHOULD  BE  ALTERED  BY  VNORM. 

FOR  EXAMPLE,  A  USER-SUPPLIED  VNORM  ROUTINE  MIGHT. . 

-SUBSTITUTE  A  MAX-NORM  OF  (V(I)*W(I))  FOR  THE  RMS-NORM,  OR 
-IGNORE  SOME  COMPONENTS  OF  V  IN  THE  NORM,  WITH  THE  EFFECT  OF 
SUPPRESSING  THE  ERROR  CONTROL  ON  THOSE  COMPONENTS  OF  Y. 


C  OTHER  ROUTINES  IN  THE  LSODES  PACKAGE. 

^  IN  ADDITION  TO  SUBROUTINE  LSODES,  THE  LSODES  PACKAGE  INCLUDES  THE  - 
C  FOLLOWING  SUBROUTINES  AND  FUNCTION  ROUTINES.. 

C,  IPREP  ACTS  AS  AN  ITERFACE  BETWEEN  LSODES  AND  PREP,  AND  ALSO  DOES 
ADJUSTING  OF  WORK  SPACE  POINTERS  AND  WORK  ARRAYS. 

PREP  IS  CALLED  BY  IPREP  TO  COMPUTE  SPARSITY  AND  DO  SPARSE  MATRIX 

C  PREPROCESSING  IF  MITER  =  1  OR  2 . 


C  JGROUP  IS  CALLED  BY  PREP  TO  COMPUTE  GROUPS  OF  JACOBIAN  COLUMN 
C  INDICES  FOR  USE  WHEN  MITER  =  2. 

C  ADJLR  ADJUSTS  THE  LENGTH  OF  REQUIRED  SPARSE  MATRIX  WORK  SPACE. 

C  IT  IS  CALLED  BY  PREP. 

C  CNTNZU  IS  CALLED  BY  PREP  AND  COUNTS  THE  NONZERO  ELEMENTS  IN  THE 

C  STRICT  UPPER  TRIANGLE  OF  J  +  J-TRANSPOSE,  WHERE  J  =  DF/DY. 

C  INTDY  COMPUTES  AN  INTERPOLATED  VALUE  OF  THE  Y  VECTOR  AT  T  =  TOUT . 

I  C  STODE  IS  THE  CORE  INTEGRATOR,  WHICH  DOES  ONE  STEP  OF  THE 

C  INTEGRATION  AND  THE  ASSOCIATED  ERROR  CONTROL. 

C  CFODE  SETS  ALL  METHOD  COEFFICIENTS  AND  TEST  CONSTANTS. 

C  PRJS  COMPUTES  AND  PREPROCESSES  THE  JACOBIAN  MATRIX  J  =  DF/DY 

C  AND  THE  NEWTON  ITERATION  MATRIX  P  =  I  -  H*L0*J. 

C  SLSS  MANAGES  SOLUTION  OF  LINEAR  SYSTEM  IN  CHORD  ITERATION. 

C  EWSET  SETS  THE  ERROR  WEIGHT  VECTOR  EWT  BEFORE  EACH  STEP. 

C  VNORM  COMPUTES  THE  WEIGHTED  R.M.S.  NORM  OF  A  VECTOR. 

C  SVCMS  AND  RSCMS  ARE  USER-CALLABLE  ROUTINES  TO  SAVE  AND  RESTORE, 

C  RESPECTIVELY,  THE  CONTENTS  OF  THE  INTERNAL  COMMON  BLOCKS. 

C  ODRV  CONSTRUCTS  A  REORDERING  OF  THE  ROWS  AND  COLUMNS  OF 

C  A  MATRIX  BY  THE  MINIMUM  DEGREE  ALGORITHM.  ODRV  IS  A 

C  DRIVER  ROUTINE  WHICH  CALLS  SUBROUTINES  MD,  MDI ,  MDM, 

C  MDP,  MDU,  AND  SRO.  SEE  REF.  2  FOR  DETAILS.  (THE  ODRV 

C  MODULE  HAS  BEEN  MODIFIED  SINCE  REF.  2,  HOWEVER.) 

C  CDRV  PERFORMS  REORDERING,  SYMBOLIC  FACTORIZATION,  NUMERICAL 

C  FACTORIZATION,  OR  LINEAR  SYSTEM  SOLUTION  OPERATIONS, 

C  DEPENDING  ON  A  PATH  ARGUMENT  I PATH.  CDRV  IS  A 

C  DRIVER  ROUTINE  WHICH  CALLS  SUBROUTINES  NROC,  NSFC, 

C  NNFC,  NNSC,  AND  NNTC.  SEE  REF.  3  FOR  DETAILS. 

C  LSODES  USES  CDRV  TO  SOLVE  LINEAR  SYSTEMS  IN  WHICH  THE 

C  COEFFICIENT  MATRIX  IS  P  =  I  -  CON*J,  WHERE  I  IS  THE 

C  IDENTITY,  CON  IS  A  SCALAR,  AND  J  IS  AN  APPROXIMATION  TO 

C  THE  JACOBIAN  DF/DY.  BECAUSE  CDRV  DEALS  WITH  ROWWISE 

C  SPARSITY  DESCRIPTIONS,  CDRV  WORKS  WITH  P-TRANSPOSE,  NOT  P. 

C  RIMACH  COMPUTES  THE  UNIT  ROUNDOFF  IN  A  MACHINE- INDEPENDENT  MANNER. 
C  XERRWV,  XSETUN,  AND  XSETF  HANDLE  THE  PRINTING  OF  ALL  ERROR 
C  MESSAGES  AND  WARNINGS.  XERRWV  IS  MACHINE - DEPENDENT . 

C  NOTE . .  VNORM  AND  RIMACH  ARE  FUNCTION  ROUTINES . 

C  ALL  THE  OTHERS  ARE  SUBROUTINES. 

C 

C  THE  INTRINSIC  AND  EXTERNAL  ROUTINES  USED  BY  LSODES  ARE.. 

C  ABS,  AMAXl,  AMINl,  FLOAT,  MAXO,  MINO,  MOD,  SIGN,  SQRT,  AND  WRITE. 

C 

C - - 

C  THE  FOLLOWING  CARD  IS  FOR  OPTIMIZED  COMPILATION  ON  LLL  COMPILERS. 
CLLL.  OPTIMIZE 

C - 

EXTERNAL  PRJS,  SLSS 

INTEGER  ILLIN,  INIT,  LYH,  LEWT,  LACOR,  LSAVF,  LWM,  LIWM, 

1  MXSTEP,  MXHNIL,  NHNIL,  NTREP,  NSLAST,  NYH,  lOWNS 
INTEGER  ICF,  lERPJ,  lERSL,  JCUR,  JSTART,  KFLAG,  L,  METH,  MITER, 

1  MAXORD,  MAXCOR,  MSBP,  MXNCF,  N,  NQ,  NST,  NFE,  NJE,  NQU 
INTEGER  IPLOST,  lESP,  ISTATC,  lYS,  IBA,  IBIAN,  IBJAN,  IBJGP, 

1  IPIAN,  IPJAN,  IPJGP,  IPIGP,  IPR,  IPC,  IPIC,  IPISP,  IPRSP,  IPA, 

2  LENYH,  LENYHM,  LENWK,  LREQ,  LRAT,  LREST,  LWMIN,  MOSS,  MSBJ, 

3  NSLJ,  NGP,  NLU,  NNZ,  NSP,  NZL,  NZU 

INTEGER  I,  II,  12,  I FLAG,  IMAX,  IMUL,  IMXER,  IPFLAG,  IPGO,  TREM; 

1  J,  KGO,  LENRAT,  LENYHT,  LENIW,  LENRW,  LENWM,  LFO ,  LIA,  LJA, 

2  LRTEM,  LWTEM,  LYHD,  LYHN,  MFl,  MORD,  MXHNLO ,  MXSTPO,  NCOLM 
REAL  TRET,  ROWNS, 

1  CCMAX,  ELO,  H,  HMIN,  HMXI,  HU,  RC,  TN,  UROUND 
REAL  CONO,  CONMIN,  CCMXJ,  PSMALL,  RBIG,  SETH 


Jo  r 


REAL  ATOLI,  AYI ,  BIG,  EWTI,  HO,  HMAX,  HMX,  RH,  RTOLI, 

1  TCRIT,  TDIST,  TNEXT,  TOL,  TOLSF,  TP,  SIZE,  SUM,  WO, 

2  RIMACH,  VNORM 
DIMENSION  MORD(2) 

LOGICAL  IHIT  . 

_C - - - 

:  THE  FOLLOWING  TWO  INTERNAL  COMMON  BLOCKS  CONTAIN 

?  (A)  VARIABLES  WHICH  ARE  LOCAL  TO  ANY  SUBROUTINE  BUT  WHOSE  VALUES  MUST 
C  BE  PRESERVED  BETWEEN  CALLS  TO  THE  ROUTINE  (OWN  VARIABLES),  AND 

-T  (B)  VARIABLES  WHICH  ARE  COMMUNICATED  BETWEEN  SUBROUTINES. 

’  THE  STRUCTURE  OF  EACH  BLOCK  IS  AS  FOLLOWS . .  ALL  REAL  VARIABLES  ARE 
C  LISTED  FIRST,  FOLLOWED  BY  ALL  INTEGERS.  WITHIN  EACH  TYPE,  THE 
X  VARIABLES  ARE  GROUPED  WITH  THOSE  LOCAL  TO  SUBROUTINE  LSODES  FIRST, 

;  THEN  THOSE  LOCAL  TO  SUBROUTINE  STODE  OR  SUBROUTINE  PRJS 
J  (NO  OTHER  ROUTINES  HAVE  OWN  VARIABLES) ,  AND  FINALLY  THOSE  USED 
C  FOR  COMMUNICATION.  THE  BLOCK  LSOOOl  IS  DECLARED  IN  SUBROUTINES 

LSODES,  I  PREP,  PREP,  IIJTDY,  STODE,  PRJS,  AND  SLSS.  THE  BLOCK  LSSOOl 
:  IS  DECLARED  IN  SUBROUTINES  LSODES,  IPREP,  PREP,  PRJS,  AND  SLSS. 

C  GROUPS  OF  VARIABLES  ARE  REPLACED  BY  DUMMY  ARRAYS  IN  THE  COMMON 
-C  DECLARATIONS  IN  ROUTINES  WHERE  THOSE  VARIABLES  ARE  NOT  USED. 


COMMON  /LSOOOl/  TRET,  ROWNS(209), 

__  1  CCMAX,  ELO,  H,  HMIN,  HMXI ,  HU,  RC,  TN,  UROUND, 

2  ILLIN,  INIT,  LYH,  LEWT,  LACOR,  LSAVF,  LWM,  LIWM, 

3  MXSTEP,  MXHNIL,  NHNIL,  NTREP,  NSLAST,  NYH,  IOWNS(6), 

4  ICF,  lERPJ,  lERSL,  JCUR,  JSTART,  KFLAG,  L,  METH,  MITER, 

—  5  MAXORD,  MAXCOR,  MSBP,  MXNCF,  N,  NQ,  NST,  NFE,  NJE,  NQU 

COMMON  /LSSOOl/  CONO ,  CONMIN,  CCMXJ,  PSMALL,  RBIG,  SETH, 

_  1  IPLOST,  lESP,  ISTATC,  lYS,  IBA,  IBIAN,  IBJAN,  IBJGP, 

2  IPIAN,  IPJAN,  IPJGP,  IPIGP,  IPR,  IPC,  IPIC,  IPISP,  IPRSP,  IPA, 

3  LENYH,  LENYHM,  LENWK,  LREQ,  LRAT,  LREST,  LWMIN,  MOSS,  MSBJ, 

4  NSLJ,  NGP,  NLU,  NN2,  NSP,  NZL,  NZU 

DATA  MORD(l) ,MORD(2) /12, 5/,  MXSTPO/500/,  MXHNLO/lO/ 

DATA  ILLIN/0/,  NTREP/O/ 

- 

IN  THE  DATA  STATEMENT  BELOW,  SET  LENRAT  EQUAL  TO  THE  RATIO  OF 
THE  WORDLENGTH  FOR  A  REAL  NUMBER  TO  THAT  FOR  AN  INTEGER.  USUALLY, 
LENRAT  =  1  FOR  SINGLE  PRECISION  AND  2  FOR  DOUBLE  PRECISION.  IF  THE 
TRUE  RATIO  IS  NOT  AN  INTEGER,  USE  THE  NEXT  SMALLER  INTEGER  ( . GE .  1) . 


DATA  LENRAT/ 1/ 

_ _ _ _ _ _ _ _ 

BLOCK  A. 

C  THIS  CODE  BLOCK  IS  EXECUTED  ON  EVERY  CALL. 

Q  IT  TESTS  ISTATE  AND  ITASK  FOR  LEGALITY  AND  BRANCHES  APPROPIATELY . 

IF  ISTATE  .GT.  1  BUT  THE  FLAG  INIT  SHOWS  THAT  INITIALIZATION  HAS 
NOT  YET  BEEN  DONE,  AN  ERROR  RETURN  OCCURS. 

C  IF  ISTATE  =  1  AND  TOUT  =  T,  JUMP  TO  BLOCK  G  AND  RETURN  IMMEDIATELY. 


IF  (ISTATE  .LT.  1  .OR.  ISTATE  .GT.  3)  GO  TO  601 
IF  (ITASK  .LT.  1  .OR.  ITASK  .GT.  5)  GO  TO  602 
IF  (ISTATE  .EQ.  1)  GO  TO  10 
IF  (INIT  .EQ.  0)  GO  TO  603 
'  IF  (ISTATE  .EQ.  2)  GO  TO  200 

GO  TO  20 
10  INIT  =  0 

IF  (TOUT  .EQ.  T)  GO  TO  430 
20  NTREP  =  0 


c - - - - - 

C  BLOCK  B- 

C  THE  NEXT  CODE  BLOCK  IS  EXECUTED  FOR  THE  INITIAL  CALL  (ISTATE  =1), 

C  OR  FOR  A  CONTINUATION  CALL  WITH  PARAMETER  CHANGES  (ISTATE  =  3) . 

C  IT  CONTAINS  CHECKING  OF  ALL  INPUTS  AND  VARIOUS  INITIALIZATIONS. 

C  IF  ISTATE  =  1,  THE  FINAL  SETTING  OF  WORK  SPACE  POINTERS,  THE  MATRIX 
C  PREPROCESSING,  AND  OTHER  INITIALIZATIONS  ARE  DONE  IN  BLOCK  C. 

C 

C  FIRST  CHECK  LEGALITY  OF  THE  NON -OPTIONAL  INPUTS  NEQ,  ITOL,  lOPT, 

C  MF,  ML,  AND  MU.  - 

C - - 

IF  (NEQ(l)  .LE.  0)  GO  TO  604 

IF  (ISTATE  .EQ.  1)  GO  TO  25 

IF  (NEQ(l)  .GT.  N)  GO  TO  605 

25  N  =  NEQ(l) 

IF  (ITOL  .LT.  1  .OR.  ITOL  .GT.  4)  GO  TO  606 

IF  (lOPT  .LT.  0  .OR.  lOPT  .GT.  1)  GO  TO  607 

MOSS  =  MF/100 

MFl  =  MF  -  100*MOSS 

METH  =  MFl/10 

MITER  =  MFl  -  10*METH 

IF  (MOSS  .LT.  0  .OR.  MOSS  .GT.  2)  GO  TO  608 

IF  (METH  .LT.  1  .OR.  METH  .GT.  2)  GO  TO  608 

IF  (MITER  .LT.  0  .OR.  MITER  .GT.  3)  GO  TO  608 

IF  (MITER  .EQ.  0  .OR.  MITER  .EQ.  3)  MOSS  =  0 

C  NEXT  PROCESS  AND  CHECK  THE  OPTIONAL  INPUTS.  - 

IF  (lOPT  .EQ.  1)  GO  TO  40 
MAXORD  =  MORD(METH) 

MXSTEP  =  MXSTPO 

MXHNIL  =  MXHNLO 

IF  (ISTATE  .EQ.  1)  HO  =  0 . OEO 

HMXI  =  O.OEO 

HMIN  =  0 . OEO 

SETH  =  0 . OEO 

GO  TO  60 

40  MAXORD  =  IWORK(5) 

IF  (MAXORD  .LT.  0)  GO  TO  611 
IF  (MAXORD  .EQ.  0)  MAXORD  =  100 
MAXORD  =  MINO (MAXORD, MORD( METH) ) 

MXSTEP  =  IWORK(6) 

IF  (MXSTEP  .LT.  0)  GO  TO  612 
IF  (MXSTEP  .EQ.  0)  MXSTEP  =  MXSTPO 
MXHNIL  =  I WORK (7) 

IF  (MXHNIL  .LT.  0)  GO  TO  613 
IF  (MXHNIL  .EQ.  0)  MXHNIL  =  MXHNLO 
IF  (ISTATE  .NE.  1)  GO  TO  50 
HO  =  RWORK(5) 

IF  ((TOUT  -  T)*H0  .LT.  O.OEO)  GO  TO  614 
50  HMAX  =  RWORK(6) 

IF  (HMAX  .LT.  O.OEO)  GO  TO  615 
HMXI  =  O.OEO 

IF  (HMAX  .GT.  O.OEO)  HMXI  =  l.OEO/HMAX 
HMIN  =  RWORK(7) 

IF  (HMIN  .LT.  O.OEO)  GO  TO  616 
SETH  =  RWORK(8) 

IF  (SETH  .LT.  O.OEO)  GO  TO  609 

C  CHECK  RTOL  AND  ATOL  FOR  LEGALITY.  - 

60  RTOLI  =  RTOL(l) 

ATOL I  =  ATOL(l) 

DO  65  I  =  1,N 


J  O  ( 


IF  (ITOL  .GE.  3)  RTOLI  =  RTOL(I) 

IF  (ITOL  .EQ.  2  .OR.  ITOL  .EQ.  4)  ATOLI  =  ATOL(I) 

—  IF  (RTOLI  .LT.  O.OEO)  GO  TO  619 

IF  (ATOLI  .LT.  O.OEO)  GO  TO  620 

65  CONTINUE 

^ - - - 

;  COMPUTE  REQUIRED  WORK  ARRAY  LENGTHS,  AS  FAR  AS  POSSIBLE,  AND  TEST 
:  THESE  AGAINST  LRW  AND  LIW.  THEN  SET  TENTATIVE  POINTERS  FOR  WORK 
C  ARRAYS.  POINTERS  TO  RWORK/IWORK  SEGMENTS.  ARE  NAMED  BY  PREFIXING  L  TO 
-r  THE  NAME  OF  THE  SEGMENT.  E.G.,  THE  SEGMENT  YH  STARTS  AT  RWORK(LYH) . 

'  SEGMENTS  OF  RWORK  (IN  ORDER)  ARE  DENOTED  WM,  YH,  SAVE,  EWT,  ACOR. 

C  IF  MITER  =  1  OR  2,  THE  REQUIRED  LENGTH  OF  THE  MATRIX  WORK  SPACE  WM 
j:  is  not  yet  known,  and  so  a  crude  minimum  VALUE  IS  USED  FOR  THE 
INITIAL  TESTS  OF  LRW  AND  LIW,  AND  YH  IS  TEMPORARILY  STORED  AS  FAR 
TO  THE  RIGHT  IN  RWORK  AS  POSSIBLE,  TO  LEAVE  THE  MAXIMUM  AMOUNT 
OF  SPACE  FOR  WM  FOR  MATRIX  PREPROCESSING.  THUS  IF  MITER  =  1  OR  2 
AND  MOSS  .NE.  2,  SOME  OF  THE  SEGMENTS  OF  RWORK  ARE  TEMPORARILY 
OMITTED,  AS  THEY  ARE  NOT  NEEDED  IN  THE  PREPROCESSING.  THESE 
C  OMITTED  SEGMENTS  ARE..  ACOR  IF  ISTATE  =  1,  EWT  AND  ACOR  IF  ISTATE  =  3 
-c  AND  MOSS  =  1,  AND  SAVE,  EWT,  AND  ACOR  IF  ISTATE  =  3  AND  MOSS  =  0. 


LRAT  =  LENRAT 

_  IF  (ISTATE  .EQ.  1)  NYH  =  N 

LWMIN  =  0 

IF  (MITER  .EQ.  1)  LWMIN  =  4*N  +  10*N/LRAT 
IF  (MITER  .EQ.  2)  LWMIN  =  4*N  +  11*N/LRAT 

-  IF  (MITER  .EQ.  3)  LWMIN  =  N  +  2 
LENYH  =  (MAXORD+1) *NYH 

LREST  =  LENYH  +  3*N 
__  LENRW  =  20  +  LWMIN  +  LREST 

IWORK(17)  =  LENRW 
LENIW  =30 

IF  (MOSS  .EQ.  0  .AND.  MITER  .NE.  0  .AND.  MITER  .NE.  3) 

~  1  LENIW  =  LENIW  +  N  +  1 

IWORK(18)  =  LENIW 
IF  (LENRW  .GT.  LRW)  GO  TO  617 
~  IF  (LENIW  .GT.  LIW)  GO  TO  618 

LIA  =  31 

IF  (MOSS  .EQ.  0  .AND.  MITER  .NE.  0  .AND.  MITER  .NE.  3) 

_  1  LENIW  =  LENIW  +  IWORK(LIA+N)  -  1 

IWORK(18)  =  LENIW 
IF  (LENIW  .GT.  LIW)  GO  TO  618 
LJA  =  LIA  +  N  +  1 

-  LIA  =  MINO (LIA,LIW) 

LJA  =  MINO  (LJA,  LIW) 

LWM  =  21 

_  IF  (ISTATE  .EQ.  1)  NQ  =  1 

NCOLM  =  MINO (NQ+l,MAXORD+2) 

LENYHM  =  NCOLM*NYH 
LENYHT  =  LENYH 

IF  (MITER  .EQ.  1  .OR.  MITER  . EQ .  2)  LENYHT  =  LENYHM 
IMUL  =  2 

IF  (ISTATE  .EQ.  3)  IMUL  =  MOSS 

-  IF  (MOSS  .EQ.  2)  IMUL  =  3 
LRTEM  =  LENYHT  +  IMUL*N 

^  LWTEM  =  LWMIN 

_  IF  (MITER  .EQ.  1  .OR.  MITER  .EQ.  2)  LWTEM  =  LRW  -  20  -  LRTEM 

LENWK  =  LWTEM 
LYHN  =  LWM  +  LWTEM 
LSAVF  =  LYHN  +  LENYHT 


LEWT  =  LSAVF  +  N 

LACOR  =  LEWT  +  N 

ISTATC  =  I STATE 

IF  (ISTATE  .EQ.  1)  GO  TO  100 

C - 

C  ISTATE  =  3.  MOVE  YH  TO  ITS  NEW  LOCATION. 

C  NOTE  THAT  ONLY  THE  PART  OF  YH  NEEDED  FOR  THE  NEXT  STEP,  NAMELY 
C  MIN(NQ+l,MAXORD+2)  COLUMNS,  IS  ACTUALLY  MOVED. 

C  A  TEMPORARY  ERROR  WEIGHT  ARRAY  EWT  IS  LOADED  IF  MOSS  =  2. 

C  SPARSE  MATRIX  PROCESSING  IS  DONE-  IN  IPREP/PREP  IF  MITER  =  1  OR  2 . 

C  IF  MAXORD  WAS  REDUCED  BELOW  NQ,  THEN  THE  POINTERS  ARE  FINALLY  SET 
C  SO  THAT  SAVE  IS  IDENTICAL  TO  YH ( * , MAXORD+2 ) . 

C - - - 

LYHD  =  LYH  -  LYHN 

IMAX  =  LYHN  -  1  +  LENYHM 

C  MOVE  YH.  BRANCH  FOR  MOVE  RIGHT,  NO  MOVE,  OR  MOVE  LEFT.  - 

IF  (LYHD)  70,80,74 
70  DO  72  I  =  LYHN, IMAX 

J  =  IMAX  +  LYHN  -  I 

72  RWORK(J)  =  RWORK(J+LYHD) 

GO  TO  80 

74  DO  76  I  =  LYHN, IMAX 
76  RWORK(I)  =  RWORK(I+LYHD) 

80  LYH  =  LYHN 

IWORK(22)  =  LYH 

IF  (MITER  .EQ.  0  .OR.  MITER  .EQ.  3)  GO  TO  92 
IF  (MOSS  .NE.  2)  GO  TO  85 

C  TEMPORARILY  LOAD  EWT  IF  MITER  =  1  OR  2  AND  MOSS  =  2.  . . 

CALL  EWSET  (N,  ITOL,  RTOL,  ATOL,  RWORK(LYH),  RWORK(LEWT)) 

DO  82  I  =  1,N 

IF  (RWORK(I+LEWT-l)  . LE .  O.OEO)  GO  TO  621 
,  82  RWORK(I+LEWT-l)  =  1 . OEO/RWORK ( I+LEWT-1) 

85  CONTINUE 

C  IPREP  AND  PREP  DO  SPARSE  MATRIX  PREPROCESSING  IF  MITER  =  1  OR  2.  . 

LSAVF  =  MINO (LSAVF, LRW) 

LEWT  =  MINO (LEWT, LRW) 

LACOR  =  MINO (LACOR, LRW) 

CALL  IPREP  (NEQ,  Y,  RWORK,  IWORK(LIA),  IWORK(LJA),  IPFLAG,  F,  JAC) 
LENRW  =  LWM  -  1  +  LENWK  +  LREST 
IWORK(17)  =  LENRW 

IF  (IPFLAG  .NE.  -1)  IWORK(23)  =  IPIAN 
IF  (IPFLAG  .NE.  -1)  IWORK(24)  =  IPJAN 
IPGO  =  -IPFLAG  +  1 

GO  TO  (90,  628,  629,  630,  631,  632,  633),  IPGO 
90  IWORK(22)  =  LYH 

IF  (LENRW  .GT.  LRW)  GO  TO  617 

C  SET  FLAG  TO  SIGNAL  PARAMETER  CHANGES  TO  STODE .  - 

92  JSTART  =  -1 

IF  (N  .EQ.  NYH)  GO  TO  200 

C  NEQ  WAS  REDUCED.  ZERO  PART  OF  YH  TO  AVOID  UNDEFINED  REFERENCES.  - 

11  =  LYH  +  L*NYH 

12  =  LYH  +  (MAXORD  +  1)*NYH  -  1 
IF  (II  .GT.  12)  GO  TO  200 

DO  95  I  =  II, 12 
95  RWORK (I)  =  O.OEO 

GO  TO  200 

C - 

C  BLOCK  C. 

C  THE  NEXT  BLOCK  IS  FOR  THE  INITIAL  CALL  ONLY  (ISTATE  =  1) . 

C  IT  CONTAINS  ALL  REMAINING  INITIALIZATIONS,  THE  INITIAL  CALL  TO  F, 


:  THE  SPARSE  MATRIX  PREPROCESSING  (MITER  =  1  OR  2),  AND  THE 
C  CALCULATION  OF  THE  INITIAL  STEP  SIZE. 

THE  ERROR  WEIGHTS  IN  EWT  ARE  INVERTED  AFTER  BEING  LOADED. 


100  CONTINUE 
_  LYH  =  LYHN 

IWORK(22)  =  LYH 
TN  =  T 
NST  =  0 

“  H  =  l.OEO 

NNZ  =  0 
NGP  =  0 
_  NZL  =  0 

NZU  =  0 

^  LOAD  THE  INITIAL  VALUE  VECTOR  IN  YH.  - 

DO  105  I  =  1,N 

■"105  RWORK(I+LYH-l)  =  Y(I) 

INITIAL  CALL  TO  F.  (LFO  POINTS  TO  YH{*,2).)  - 

LFO  =  LYH  +  NYH 

—  CALL  F  (NEQ,  T,  Y,  RWORK(LFO)) 

NFE  =  1 

LOAD  AND  INVERT  THE  EWT  ARRAY.  (H  IS  TEMPORARILY  SET  TO  1.0.)  - 

_  CALL  EWSET  (N,  ITOL,  RTOL,  ATOL,  RWORK{LYH),  RWORK(LEWT)) 

DO  110  I  =  1,N 

IF  (RWORK(I+LEWT-l)  . LE .  O.OEO)  GO  TO  621 
110  RWORK(I+LEWT-l)  =  1 . OEO/RWORK (I+LEWT-1) 

—  IF  (MITER  .EQ.  0  .OR.  MITER  . EQ .  3)  GO  TO  120 

IPREP  AND  PREP  DO  SPARSE  MATRIX  PREPROCESSING  IF  MITER  =  1  OR  2 .  - 

LACOR  =  MINO (LACOR,LRW) 

~  CALL  IPREP  (NEQ,  Y,  RWORK,  IWORK(LIA),  IWORK(LJA),  IPFLAG,  F,  JAC) 

LENRW  =  LWM  -  1  +  LENWK  +  LREST 
IWORK(17)  =  LENRW 

_  IF  (IPFLAG  .NE.  -1)  IWORK(23)  =  IPIAN 

IF  (IPFLAG  .NE.  -1)  IWORK(24)  =  IPJAN 
IPGO  =  -IPFLAG  +  1 

GO  TO  (115,  628,  629,  630,  631,  632,  633),  IPGO 
—115  IWORK(22)  =  LYH 

IF  (LENRW  .GT.  LRW)  GO  TO  617 

u  CHECK  TCRIT  FOR  LEGALITY  (ITASK  =  4  OR  5).  - 

_120  CONTINUE 

IF  (ITASK  .NE.  4  .AND.  ITASK  .NE.  5)  GO  TO  125 
TCRIT  =  RWORK (1) 

IF  ((TCRIT  -  TOUT)* (TOUT  -  T)  .LT.  O.OEO)  GO  TO  625 
~  IF  (HO  .NE.  O.OEO  .AND.  (T  +  HO  -  TCRIT) *H0  . GT .  O.OEO) 

1  HO  =  TCRIT  -  T 

C  INITIALIZE  ALL  REMAINING  PARAMETERS.  - : - 

~12S  UROUND  =  R1MACH(4) 

JSTART  =  0 

IF  (MITER  .NE.  0)  RWORK (LWM)  =  SQRT (UROUND) 

_  MSBJ  =  50 

NSLJ  =  0 
CCMXJ  =  0.2E0 
PSMALL  =  1000 . 0E0*UROUND 

—  RBIG  =  0  .  OlEO/PSMTUiL 

^  NHNIL  =0  ‘  ' 

NJE  =  0 
_  NLU  =  0 

NSLAST  =  0 
HU  =  O.OEO 


CCMAX  =  0.3E0 
MAXCOR  =  3 
MSBP  =  20 
MXNCF  =  10 

C - 

C  THE  CODING  BELOW  COMPUTES  THE  STEP  SIZE,  HO,  TO  BE  ATTEMPTED  ON  THE 
C  FIRST  STEP,  UNLESS  THE  USER  HAS  SUPPLIED  A  VALUE  FOR  THIS. 

C  FIRST  CHECK  THAT  TOUT  -  T  DIFFERS  SIGNIFICANTLY  FROM  ZERO. 

C  A  SCALAR  TOLERANCE  QUANTITY  TOL  IS  COMPUTED,  AS  MAX(RTOL(I)) 

C  IF  THIS  IS  POSITIVE,  OR  MAX  (ATOL"(  I ) /ABS  (Y  ( I )  )  )  OTHERWISE,  ADJUSTED 
C  SO  AS  TO  BE  BETWEEN  100*UROUND  AND  l.OE-3. 

C  THEN  THE  COMPUTED  VALUE  HO  IS  GIVEN  BY. . 

C  NEQ 

C  H0**2  =  TOL  /  (  W0**-2  +  (1/NEQ)  *  SUM  {  F ( I ) /YWT { I )  ) **2  ) 

C  1 

C  WHERE  WO  =  MAX  (  ABS(T),  ABS (TOUT)  ), 

C  F(I)  =  I-TH  COMPONENT  OF  INITIAL  VALUE  OF  F, 

C  YWT (I)  =  EWT(I)/TOL  (A  WEIGHT  FOR  Y(I)). 

C  THE  SIGN  OF  HO  IS  INFERRED  FROM  THE  INITIAL  VALUES  OF  TOUT  AND  T. 

C - 

LFO  =  LYH  +  NYH 

IF  (HO  .NE.  O.OEO)  GO  TO  180 

TDIST  =  ABS (TOUT  -  T) 

WO  =  AMAXl (ABS (T) , ABS (TOUT) ) 

IF  (TDIST  .LT.  2 . OEO*UROUND*WO)  GO  TO  622 
TOL  =  RTOL(l) 

IF  (ITOL  .LE.  2)  GO  TO  140 
DO  130  I  =  1,N 

130  TOL  =  AMAXl (TOL,RTOL(I) ) 

140  IF  (TOL  .GT.  O.OEO)  GO  TO  160 
ATOLI  =  ATOL(l) 

DO  150  I  =  1,N 

IF  (ITOL  .EQ.  2  .OR.  ITOL  . EQ .  4)  ATOLI  =  ATOL(I) 

AYI  =  ABS (Y (I) ) 

IF  (AYI  .NE.  O.OEO)  TOL  =  AMAXl (TOL, ATOLI /AYI ) 

150  CONTINUE 

160  TOL  =  AMAXl (TOL, 100 . 0E0*UROUND) 

TOL  =  AMINl (TOL, 0 . OOIEO) 

SUM  =  VNORM  (N,  RWORK(LFO),  RWORK(LEWT)) 

SUM  =  1 . OEO/ (TOL*W0*W0)  +  TOL*SUM**2 
HO  =  1.0E0/SQRT(SUM) 

HO  =  AMINl (HO, TDIST) 

HO  =  SIGN(H0,TOUT-T) 

C  ADJUST  HO  IF  NECESSARY  TO  MEET  HMAX  BOUND.  - 

180  RH  =  ABS(H0)*HMXI 

IF  (RH  .GT.  l.OEO)  HO  =  HO/RH 

C  LOAD  H  WITH  HO  AND  SCALE  YH(*,2)  BY  HO.  - 

H  =  HO 

DO  190  I  =  1,N 

190  RWORK(I+LF0-l)  =  H0*RWORK ( I+LFO-1) 

GO  TO  270 

I  c - 

^  C  BLOCK  D. 

C  THE  NEXT  CODE  BLOCK  IS  FOR  CONTINUATION  CALLS  ONLY  (ISTATE  =  2  OR  3 ) 
C  AND  IS  TO  CHECK  STOP  CONDITIONS  BEFORE  TAKING  A  STEP. 

C---: - - - 

200  NSLAST  =  NST 

GO  TO  (210,  250,  220,  230,  240),  ITASK 
210  IF  ( (TN  -  TOUT) *H  . LT .  O.OEO)  GO  TO  250 

CALL  INTDY  (TOUT,  0,  RWORK(LYH),  NYH,  Y,  IFLAG) 


1  O  (  )">  c 


IF  (IFLAG  .NE.  0)  GO  TO  627 
T  =  .TOUT 
GO  TO  420 

220  ,TP  =  TN  -  HU*{1.0E0  +  100 . OEO *UROUND) 

IF  ((TP  -  TOUT)*H  .GT.  O.OEO)  GO  TO  623 
IF  ( (TN  -  TOUT)*H  .LT.  O.OEO)  GO  TO  250 
GO  TO  400 

'230  TCRIT  =  RWORK(l) 

IF  ( (TN  -  TCRIT) *H  . GT .  O.OEO)  GO  TO  624 
]  IF  ((TCRIT  -  TOUT) *H  .LT.  O.UEO)  GO  TO  625 

IF  ( (TN  -  TOUT)*H  .LT.  O.OEO)  GO  TO  245 
CALL  INTDY  (TOUT,  0,  RWORK(LYH),  NYH,  Y,  IFLAG) 
,  IF  (IFLAG  .NE.  0)  GO  TO  627 

T  =  TOUT 
GO  TO  420 

240  TCRIT  =  RWORK(l) 

IF  ( (TN  -  TCRIT) *H  .GT.  O.OEO)  GO  TO,  624 
245  HMX  =  ABS(TN)  +  ABS (H) 

IHIT  =  ABS(TN  -  TCRIT)  . LE .  100 . 0E0*UROUND*HMX 
IF  (IHIT)  GO  TO  400 

TNEXT  =  TN  +  H*(1.0E0  +  4 . 0E0*UROUND) 

IF  ((TNEXT  -  TCRIT) *H  . LE .  O.OEO)  GO  TO  250 
H  =  (TCRIT  -  TN)*(1.0E0  -  4 . 0E0*UROUND) 

'  IF  (ISTATE  .EQ.  2)  JSTART  =  -2 


C  BLOCK  E. 

^  THE  NEXT  BLOCK  IS  NORMALLY  EXECUTED  FOR  ALL  CALLS  AND  CONTAINS 
THE  CALL  TO  THE  ONE -STEP  CORE  INTEGRATOR  STODE. 

THIS  IS  A  LOOPING  POINT  FOR  THE  INTEGRATION  STEPS. 

FIRST  CHECK  FOR  TOO  MANY  STEPS  BEING  TAKEN,  UPDATE  EWT  (IF  NOT  AT 
START  OF  PROBLEM) ,  CHECK  FOR  TOO  MUCH  ACCURACY  BEING  REQUESTED,  AND 
CHECK  FOR  H  BELOW  THE  ROUNDOFF  LEVEL  IN  T. 


250  CONTINUE 

—  IF  ( (NST-NSLAST)  . GE .  MXSTEP)  GO  TO  500 

CALL  EWSET  (N,  ITOL,  RTOL,  ATOL,  RWORK(LYH),  RWORK(LEWT)) 
DO  260  I  =  1,N 

__  IF  (RWORK(I+LEWT-l)  .  LE .  O.OEO)  GO  TO  510 

260  RWORK(I+LEWT-l)  =  1 . OEO/RWORK (I+LEWT-1) 

270  TOLSF  =  UROUND*VNORM  (N,  RWORK(LYH),  RWORK(LEWT)) 

IF  (TOLSF  .LE.  l.OEO)  GO  TO  280 

—  TOLSF  =  TOLSF*2.0E0 

"1  IF  (NST  .EQ.  0)  GO  TO  626 

GO  TO  520 

_280  IF  ( (TN  +  H)  .NE.  TN)  GO  TO  290 
NHNIL  =  NHNIL  +  1 
IF  (NHNIL  .GT.  MXHNIL)  GO  TO  290 


CALL  XERRWV(50HLSODES 

WARNING. . INTERNAL  T 

(=R1) 

AND 

H  (=R2)  ARE, 

1  50,  101, 

1,  0,  0, 

0, 

0,  O.OEO,  O.OEO) 

CALL  XERRWV 
1  60H 

( 

SUCH  THAT 

IN 

THE  MACHINE,  T  +  H  = 

T  ON 

THE 

NEXT 

STEP  , 

— 

1  60,  101, 

1,  0,  0, 

0, 

0,  O.OEO,  O.OEO) 

> 

CALL  XERRWV (5 OH 

(H 

=  STEP  SIZE) .  SOLVER 

WILL 

CONTINUE 

^ANYWAY, 

1  50,  101, 

1,  0,  0, 

0, 

2,  TN,  H) 

IF  (NHNIL  .LT.  MXHNIL)  GO  TO  290 

CALL  XERRWV(50HLSODES--  ABOVE  WARNING  HAS  BEEN  ISSUED  II  TIMES. 

1  50,  102,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

CALL  XERRWV(50H  IT  WILL  NOT  BE  ISSUED  AGAIN  FOR  THIS  PROBLEM, 


ooonn  ooo 


1  50,  102,  1,  1,  MXHNIL,  0,  0,  O.OEO,  O.OEO) 

290  CONTINUE 


CALL  STODE ( NEQ , Y , YH , NYH , YH , EWT , SAVE , ACOR , WM , WM , F , JAC , PR JS , SLSS ) 


CALL  STODE  (NEQ,  Y,  RWORK(LYH),  NYH,  RWORK(LYH),  RWORK{LEWT), 

1  RWORK{LSAVF) ,  RWORK ( LACOR ) ,  RWORK{LWM),  RWORK(LWM), 

2  F,  JAC,  PRJS,  SLSS) 

KGO  =  1  -  KFLAG 

GO  TO  (300,  530,  540),  KGO  " 


BLOCK  F. 

THE  FOLLOWING  BLOCK  HANDLES  THE  CASE  OF  A  SUCCESSFUL  RETURN  FROM  THE 
CORE  INTEGRATOR  (KFLAG  =  0) .  TEST  FOR  STOP  CONDITIONS. 


300  INIT  =  1 

GO  TO  (310,  400,  330,  340,  350),  ITASK 

C  ITASK  =1.  IF  TOUT  HAS  BEEN  REACHED,  INTERPOLATE.  - 

310  IF  ( (TN  -  TOUT) *H  .LT.  O.OEO)  GO  TO  250 

CALL  INTDY  (TOUT,  0,  RWORK (LYH),  NYH,  Y,  I FLAG) 

T  =  TOUT 
GO  TO  420 

C  ITASK  =  3.  JUMP  TO  EXIT  IF  TOUT  WAS  REACHED.  - 

330  IF  ( (TN  -  TOUT)*H  . GE .  O.OEO)  GO  TO  400 
GO  TO  250 

C  ITASK  =  4.  SEE  IF  TOUT  OR  TCRIT  WAS  REACHED.  ADJUST  H  IF  NECESSARY. 
340  IF  ( (TN  -  TOUT) *H  .LT.  O.OEO)  GO  TO  345 

CALL  INTDY  (TOUT,  0,  RWORK (LYH),  NYH,  Y,  IFLAG) 

T  =  TOUT 
GO  TO  420 

345  HMX  =  ABS(TN)  +  ABS (H) 

IHIT  =  ABS(TN  -  TCRIT)  .LE.  100 . 0E0*UROUND*HMX 
IF  (IHIT)  GO  TO  400 

TNEXT  =  TN  +  H*(1.0E0  +  4 . OEO *UROUND) 

IF  ((TNEXT  -  TCRIT) *H  . LE .  O.OEO)  GO  TO  250 
H  =  (TCRIT  -  TN)*(1.0E0  -  4 . OEO *UROUND) 

JSTART  =  -2 
GO  TO  250 

C  ITASK  =  5.  SEE  IF  TCRIT  WAS  REACHED  AND  JUMP  TO  EXIT.  - 

350  HMX  =  ABS(TN)  +  ABS (H) 

IHIT  =  ABS(TN  -  TCRIT)  . LE .  100 . 0E0*UROUND*HMX 

C - 

C  BLOCK  G. 

C  THE  FOLLOWING  BLOCK  HANDLES  ALL  SUCCESSFUL  RETURNS  FROM  LSODES . 

C  IF  ITASK  .NE.  1,  Y  IS  LOADED  FROM  YH  AND  T  IS  SET  ACCORDINGLY. 

C  ISTATE  IS  SET  TO  2,  THE  ILLEGAL  INPUT  COUNTER  IS  ZEROED,  AND  THE 
C  OPTIONAL  OUTPUTS  ARE  LOADED  INTO  THE  WORK  ARRAYS  BEFORE  RETURNING. 

C  IF  ISTATE  =  1  AND  TOUT  =  T,  THERE  IS  A  RETURN  WITH  NO  ACTION  TAKEN, 

C  EXCEPT  THAT  IF  THIS  HAS  HAPPENED  REPEATEDLY,  THE  RUN  IS  TERMINATED. 

C - - - 

400  DO  410  I  =  1,N 

410  Y(I)  =  RWORK(I+LYH-l) 

T  =  TN 

IF  (ITASK  .NE.  4  .AND.  ITASK  .NE.  5)  GO  TO  420 
IF  (IHIT)  T  =  TCRIT 
420  ISTATE  =2 
ILLIN  =  0 
RWORK (11)  =  HU 
RWORK (12)  =  H 
RWORK (13)  =  TN 


IWORK(ll)  =  NST 
IWORK{12)  =  NFE 
IWORK(13)  =  NJE 
IWORK(14)  =  NQU 
IWORK(15)  =  NQ 

_  IWORK(19)  =  NNZ 

IWORK(20)  =  NGP 
IWORK(21)  =  NLU 
IWORK(25)  =  NZL 
IWORK(26)  =  NZU 
RETURN 

C 

_430  NTREP  =  NTREP  +1 

IF  (NTREP  .LT.  5)  RETURN 
CALL  XERRWV( 

_  1  60HLSODES--  REPEATED  CALLS  WITH  ISTATE  =  1  AND  TOUT  =  T  (=R1) 

1  60,  301,  1,  0,  0,  0,  1,  T,  O.OEO) 

GO  TO  800 

C - 

BLOCK  H. 

THE  FOLLOWING  BLOCK  HANDLES  ALL  UNSUCCESSFUL  RETURNS  OTHER  THAN 

THOSE  FOR  ILLEGAL  INPUT.  FIRST  THE  ERROR  MESSAGE  ROUTINE  IS  CALLED. 

C  IF  THERE  WAS  AN  ERROR  TEST  OR  CONVERGENCE  TEST  FAILURE,  IMXER  IS  SET. 

THEN  Y  IS  LOADED  FROM  YH,  T  IS  SET  TO  TN,  AND  THE  ILLEGAL  INPUT 

COUNTER  ILLIN  IS  SET  TO  0.  THE  OPTIONAL  OUTPUTS  ARE  LOADED  INTO 

C  THE  WORK  ARRAYS  BEFORE  RETURNING. 


.  THE  MAXIMUM  NUMBER  OF  STEPS  WAS  TAKEN  BEFORE  REACHING  TOUT.  . . - 

500  CALL  XERRWV(50HLSODES--  AT  CURRENT  T  {=R1),  MXSTEP  (=11)  STEPS 
_  1  50,  201,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

CALL  XERRWV(50H  TAKEN  ON  THIS  CALL  BEFORE  REACHING  TOUT 

1  50,  201,  1,  1,  MXSTEP,  0,  1,  TN,  O.OEO) 

__  ISTATE  =  -1 

GO  TO  580 

EWT(I)  .LE.  0.0  FOR  SOME  I  (NOT  AT  START  OF  PROBLEM).  . . . . . 

510  EWTI  =  RWORK(LEWT+I-l) 

—  CALL  XERRWV(50HLSODES--  AT  T  (=R1) ,  EWT(Il)  HAS  BECOME  R2  . LE .  0., 

1  50,  202,  1,  1,  I,  0,  2,  TN,  EWTI) 

ISTATE  =  -6 
_  GO  TO  580 

TOO  MUCH  ACCURACY  REQUESTED  FOR  MACHINE  PRECISION.  - 

520  CALL  XERRWV(50HLSODES--  AT  T  (=R1),  TOO  MUCH  ACCURACY  REQUESTED  , 
1  50,  203,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

—  CALL  XERRWV(50H  FOR  PRECISION  OF  MACHINE..  SEE  TOLSF  (=R2)  , 

1  50,  203,  1,  0,  0,  0,  2,  TN,  TOLSF) 

RWORK(14)  =  TOLSF 
_  ISTATE  =  -2 

GO  TO  580 

v-  KFLAG  =  -1.  ERROR  TEST  FAILED  REPEATEDLY  OR  WITH  ABS(H)  =  HMIN.  - - 

__530  CALL  XERRWV(50HLSODES--  AT  T(=R1)  AND  STEP  SIZE  H(=R2),  THE  ERROR, 
1  50,  204,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

CALL  XERRWV(50H  TEST  FAILED  REPEATEDLY  OR  WITH  ABS (H)  =  HMIN, 

1  50,  204,  1,  0,  0,  0,  2,  TN,  H) 

—  ISTATE  =  -4 

GO  TO  560 

KFLAG  =  -2.  CONVERGENCE  FAILED  REPEATEDLY  OR  WITH  ABS(H)  =  HMIN.  - 

_340  CALL  XERRWV(50HLSODES--  AT  T  (=R1)  AND  STEP  SIZE  H  (=R2),  THE 
1  50,  205,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

CALL  XERRWV(50H  CORRECTOR  CONVERGENCE  FAILED  REPEATEDLY 

1  50,  205,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 


CALL  XERRWVOOH  OR  WITH  ABS  (H)  =  HMIN 

1  30,  205,  1,  0,  0,  0,  2,  TN,  H) 

I STATE  =  -5 

C  COMPUTE  IMXER  IF  RELEVANT.  - 

560  BIG  =  O.OEO 
IMXER  =  1 
DO  570  I  =  1,N 

SIZE  =  ABS (RWORK(I+LACOR-l) *RWORK(I+LEWT-l) ) 
IF  (BIG  .GE.  SIZE)  GO  TO  570 
BIG  =  SIZE 
IMXER  =  I 
570  CONTINUE 

IWORK(16)  =  IMXER 

C  SET  Y  VECTOR,  T,  ILLIN,  AND  OPTIONAL  OUTPUTS.  - 

580  DO  590  I  =  1,N 

590  Y{I)  =  RWORK(I+LYH-l) 

T  =  TN 


ILLIN  =  0 

RWORK(ll) 

= 

HU 

RWORK(12) 

H 

RWORK(13) 

= 

TN 

IWORK(ll) 

S 

NST 

IWORK(12) 

= 

NFE 

IWORK(13) 

= 

NJE 

IWORK ( 14 ) 

= 

NQU 

IWORK(15) 

= 

NQ 

IWORK (19) 

= 

NNZ 

IWORK (20) 

= 

NGP 

IWORK (21) 

= 

NLU 

IWORK (25) 

= 

NZL 

IWORK (26) 

ss 

NZU 

RETURN 

C  BLOCK  I . 

C  THE  FOLLOWING  BLOCK  HANDLES  ALL  ERROR  RETURNS  DUE  TO  ILLEGAL  INPUT 
C  (ISTATE  =  -3) ,  AS  DETECTED  BEFORE  CALLING  THE  CORE  INTEGRATOR. 

C  FIRST  THE  ERROR  MESSAGE  ROUTINE  IS  CALLED.  THEN  IF  THERE  HAVE  BEEN 
C  5  CONSECUTIVE  SUCH  RETURNS  JUST  BEFORE  THIS  CALL  TO  THE  SOLVER, 

C  THE  RUN  IS  HALTED. 


601  CALL  XERRWV(30HLSODES--  ISTATE  (=11)  ILLEGAL  , 

1  30,  1,  1,  1,  ISTATE,  0,  0,  O.OEO,  O.OEO) 

GO  TO  700 

602  CALL  XERRWV(30HLSODES--  ITASK  (=11)  ILLEGAL  , 

1  30,  2,  1,  1,  ITASK,  0,  0,  O.OEO,  O.OEO) 

GO  TO  700 

603  CALL  XERRWV(50HLSODES--  ISTATE  . GT .  1  BUT  LSODES  NOT  INITIALIZED  , 

1  50,  3,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

GO  TO  700 

604  CALL  XERRWV(30HLSODES--  NEQ  (=11)  .LT.  1 

1  30,  4,  1,  1,  NEQ(l),  0,  0,  O.OEO,  O.OEO) 

GO  TO  700 

605  CALL  XERRWV(50HLSODES--  ISTATE  =  3  AND  NEQ  INCREASED  (II  TO  12) 

1  50,  5,  1,  2,  N,  NEQ(l),  0,  O.OEO,  O.OEO) 

GO  TO  700 

606  CALL  XERRWV(30HLSODES--  ITOL  (=11)  ILLEGAL 

1  30,  6,  1,  1,  ITOL,  0,  0,  O.OEO,  O.OEO) 

GO  TO  700 

607  CALL  XERRWV(30HLSODES--  lOPT  (=11)  ILLEGAL 

1  30,  7,  1,  1,  lOPT,  0,  0,  O.OEO,  O.OEO) 


GO  TO  700 

608  CALL  XERRWV(30HLSODES--  MF  (=11)  ILLEGAL 

—  1  30,  8,  1,  1,  MF,  0,  0,  O.OEO,  O.OEO) 

GO  TO  700 

609  CALL  XERRWV(30HLSODES--  SETH  (=R1)  .LT.  0,0  , 

_  1  30,  9,  1,  0,  0,  0,  1,  SETH,  O.OEO) 

GO  TO  700 

611  CALL  XERRWV(30HLSODES--  MAXORD  (=11)  .LT.  0  , 

1  30,  11,  1,  1,  MAXORD,  0,  0,  O.OEO,  O.OEO) 

~  GO  TO  700 

612  CALL  XERRWV(30HLSODES--  MXSTEP  (=11)  .LT.  0  , 

1  30,  12,  1,  1,  MXSTEP,  0,  0,  O.OEO,  O.OEO) 

_  GO  TO  700 

'613  CALL  XERRWV(3  0HLSODES--  MXHNIL  (  =  11)  .LT.  0  , 

1  30,  13,  1,  1,  MXHNIL,  0,  0,  O.OEO,  O.OEO) 

GO  TO  700 

”614  CALL  XERRWV(40HLSODES--  TOUT  (=R1)  BEHIND  T  (=R2) 

1  40,  14,  1,  0,  0,  0,  2,  TOUT,  T) 

CALL  XERRWV(50H  INTEGRATION  DIRECTION  IS  GIVEN  BY  HO  (=R1) 

—  1  50,  14,  1,  0,  0,  0,  1,  HO,  O.OEO) 

GO  TO  700 

615  CALL  XERRWV(30HLSODES--  HMAX  (=R1)  . LT .  0.0  , 

_  1  30,  15,  1,  0,  0,  0,  1,  HMAX,  O.OEO) 

GO  TO  700 

616  CALL  XERRWV(30HLSODES--  HMIN  (=R1)  .LT.  0.0  , 

1  30,  16,  1,  0,  0,  0,  1,  HMIN,  O.OEO) 

GO  TO  700 

617  CALL  XERRWV(50HLSODES--  RWORK  LENGTH  IS  INSUFFICIENT  TO  PROCEED.  , 

1  50,  17,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

—  CALL  XERRWV ( 

1  60H  LENGTH  NEEDED  IS  .GE.  LENRW  (=11),  EXCEEDS  LRW  (=12), 

1  60,  17,  1,  2,  LENRW,  LRW,  0,  O.OEO,  O.OEO) 

_  GO  TO  700 

618  CALL  XERRWV ( 5 OHLSODES--  IWORK  LENGTH  IS  INSUFFICIENT  TO  PROCEED.  , 

y  1  50,  18,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

CALL  XERRWV ( 

—  1  60H  LENGTH  NEEDED  IS  . GE .  LENIW  (=11),  EXCEEDS  LIW  (=12), 

1  60,  18,  1,  2,  LENIW,  LIW,  0,  O.OEO,  O.OEO) 

GO  TO  700 

_619  CALL  XERRWV ( 4 OHLSODES--  RTOL(Il)  IS  R1  .LT.  0.0 
1  40,  19,  1,  1,  I,  0,  1,  RTOLI,  O.OEO) 

GO  TO  700 

620  CALL  XERRWV ( 4 OHLSODES--  ATOL(Il)  IS  R1  .LT.  0.0 

1  40,  20,  1,  1,  I,  0,  1,  ATOLI,  O.OEO) 

GO  TO  700 

621  EWTI  =  RWORK (LEWT+I-1) 

—  CALL  XERRWV ( 4 OHLSODES--  EWT(Il)  IS  R1  .LE.  0.0 

1  40,  21,  1,  1,  I,  0,  1,  EWTI,  O.OEO) 

GO  TO  700 
_622  CALL  XERRWV ( 

1  60HLSODES--  TOUT  (=R1)  TOO  CLOSE  TO  T(=R2)  TO  START  INTEGRATION, 
1  60,  22,  1,  0,  0,  0,  2,  TOUT,  T) 

GO  TO  700 
-^23  CALL  XERRWV ( 

j  .  1  60HLSODES--  ITASK  =11  AND  TOUT  (=R1)  BEHIND  TCUR  -  HU  (='  R2)  -  , 

^  1  60,  23,  1,  1,  ITASK,  0,  2,  TOUT,  TP) 

_  GO  TO  700 
324  CALL  XERRWV ( 

1  60HLSODES--  ITASK  =  4  OR  5  AND  TCRIT  (=R1)  BEHIND  TCUR  (=R2) 

1  60,  24,  1,  0,  0,  0,  2,  TCRIT,  TN) 


GO  TO  700 

625  CALL  XERRWV( 

1  60HLSODES--  ITASK  =  4  OR  5  AND  TCRIT  (=R1)  BEHIND  TOUT  (=R2) 

1  60,  25  v  1,  0,  0,  0,  2,  TCRIT,  TOUT) 

GO  TO  700 

626  CALL  XERRWV(50HLSODES--  AT  START  OF  PROBLEM,  TOO  MUCH  ACCURACY 

1  50,  26,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

CALL  XERRWV( 

1  6 OH  REQUESTED  FOR  PRECISION  OF  MACHINE..  SEE  TOLSF  (=R1)  , 

1  60,  26,  1,  0,  0,  0,  1,  TOLSF,  O.OEO) 

RWORK ( 14 )  =  TOLSF 
GO  TO  700 

627  CALL  XERRWV(50HLSODES--  TROUBLE  FROM  INTDY.  ITASK  =  II,  TOUT  =  Rl, 

1  50,  27,  1,  1,  ITASK,  0,  1,  TOUT,  O.OEO) 

GO  TO  700 

628  CALL  XERRWV( 

1  60HLSODES--  RWORK  LENGTH  INSUFFICIENT  (FOR  SUBROUTINE  PREP) . 

1  60,  28,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

CALL  XERRWV( 

1  60H  LENGTH  NEEDED  IS  .GE.  LENRW  (=11),  EXCEEDS  LRW  (=12), 

1  60,  28,  1,  2,  LENRW,  LRW,  0,  O.OEO,  O.OEO) 

GO  TO  700 

629  CALL  XERRWV( 

1  60HLSODES--  RWORK  LENGTH  INSUFFICIENT  (FOR  SUBROUTINE  JGROUP) .  , 

1  60,  29,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

CALL  XERRWV( 

1  60H  LENGTH  NEEDED  IS  . GE .  LENRW  (=11),  EXCEEDS  LRW  (=12), 

1  60,  29,  1,  2,  LENRW,  LRW,  0,  O.OEO,  O.OEO) 

GO  TO  700 

630  CALL  XERRWV( 

1  60HLSODES--  RWORK  LENGTH  INSUFFICIENT  (FOR  SUBROUTINE  ODRV) . 

1  60,  30,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

CALL  XERRWV( 

1  6 OH  LENGTH  NEEDED  IS  .GE.  LENRW  (=11),  EXCEEDS  LRW  (=12), 

1  60,  30,  1,  2,  LENRW,  LRW,  0,  O.OEO,  O.OEO) 

GO  TO  700 

631  CALL  XERRWV( 

1  60HLSODES--  ERROR  FROM  ODRV  IN  YALE  SPARSE  MATRIX  PACKAGE 
1  60,  31,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

IMUL  =  (lYS  -  1)/N 

IREM  =  lYS  -  IMUL*N 
CALL  XERRWV( 

1  60H  AT  T  (=R1) ,  ODRV  RETURNED  ERROR  FLAG  =  I1*NEQ  +12. 

1  60,  31,  1,  2,  IMUL,  IREM,  1,  TN,  O.OEO) 

GO  TO  700 

632  CALL  XERRWV( 

1  60HLSODES--  RWORK  LENGTH  INSUFFICIENT  (FOR  SUBROUTINE  CDRV) . 

1  60,  32,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

CALL  XERRWV( 

1  6 OH  LENGTH  NEEDED  IS  .GE.  LENRW  (=11),  EXCEEDS  LRW  (=12), 

1  60,  32,  1,  2,  LENRW,  LRW,  0,  O.OEO,  O.OEO) 

GO  TO  700 

633  CALL  XERRWV( 

1  60HLSODES--  ERROR  FROM  CDRV  IN  YALE  SPARSE  MATRIX  PACKAGE 
1  60,  33,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

IMUL  =  ( lYS  -  1 ) /N 
IREM  =  lYS  -  IMUL*N 
CALL  XERRWV( 

1  6 OH  AT  T  (=R1) ,  CDRV  RETURNED  ERROR  FLAG  =  I1*NEQ  +12. 

1  60,  33,  1,  2,  IMUL,  IREM,  1,  TN,  O.OEO) 


DESCRIPTORS 


IF  (IMUL  .EQ,  2)  CALL  XERRWV ( 

1  60H  DUPLICATE  ENTRY  IN  SPARSITY  STRUCTURE 

1  60,  33,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

IF  (IMUL  .EQ.  3  .OR.  IMUL  .EQ.  6)  CALL  XERRWV ( 

1  6 OH  INSUFFICIENT  STORAGE  FOR  NSFC  (CALLED  BY  CDRV) 

1  60,  33,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

700  IF  (ILLIN  .EQ.  5)  GO  TO  710 
ILLIN  =  ILLIN  +  1 

:  ISTATE  =  -3  - 

RETURN 

710  CALL  XERRWV ( 5 OHLSODES--  REPEATED  OCCURRENCES  OF  ILLEGAL  INPUT 
1  50,  302,  1,  0,  0,  0,  0,  O.OEO,  O.OEO) 

800  CALL  XERRWV ( 5 OHLSODES--  RUN  ABORTED..  APPARENT  INFINITE  LOOP 
1  50,  303,  2,  0,  0,  0,  0,  O.OEO,  O.OEO) 

RETURN 


END 


END  OF  SUBROUTINE  LSODES 


